{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read timeseries HDF5 file and save data into text file for Comsol ... ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Go to directory: /Users/yunjunz/insarlab/Galapagos/AlcedoSenDT128/mintpy/geo\n"
     ]
    }
   ],
   "source": [
    "# read / save displacement data from time-series file\n",
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mintpy.utils import readfile\n",
    "\n",
    "work_dir = os.path.expanduser('~/insarlab/Galapagos/AlcedoSenDT128/mintpy/geo')\n",
    "os.chdir(work_dir)\n",
    "print('Go to directory:', work_dir)\n",
    "\n",
    "date12 = '20160101_20171022'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 20160101_20171022 from file geo_timeseries_ECMWF_ramp_demErr.h5\n",
      "converting range to phase\n",
      "write 20160101_20171022.unw\n",
      "write 20160101_20171022.unw.rsc\n",
      "masking 20160101_20171022 from 20160101_20171022.unw ...\n",
      "write 20160101_20171022_msk.unw\n",
      "write 20160101_20171022_msk.unw.rsc\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "# 1. save interferogram file in ROI_PAC format from time-series file\n",
    "!save_roipac.py geo_timeseries_ECMWF_ramp_demErr.h5 {date12} -o {date12+'.unw'}\n",
    "!mask.py {date12+'.unw'} -m geo_maskTempCoh.h5 -o {date12+'_msk.unw'}\n",
    "#!subset.py {date12+'_msk.unw'} --lat -1.1 -0.7 --lon -91.55 -91.0 -o {date12+'_msk_sub.unw'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3cAAADSCAYAAADpCjnWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXmcHFd177/nVvUyM9otL/ImWV4h5pmAQxaSz0vYAgRikkDY8tgxPHAAW7axDTYG22BsLYEQMAYCTgiPLZtJHAgESCBhs4EEbONNC5YtbGvXbN1ddc/7497qrh6NNCNpejQzOt/Ppz5dfetW1e2emTP3V+fcc0RVMQzDMAzDMAzDMGY37nAPwDAMwzAMwzAMwzh0TNwZhmEYhmEYhmHMAUzcGYZhGIZhGIZhzAFM3BmGYRiGYRiGYcwBTNwZhmEYhmEYhmHMAUzcGYZhGIZhGIZhzAFM3B0GROQmEbnyEK/x2yKyearGZBiGMZ2IyKCIrJzC66mInDZV1zMMwygQkX8RkVdO4fW+KSKvm6rrGUYZE3c9QEQ2isiIiOwRkZ0i8l8i8kYRcQCq+kZVveZwj3MmYYbOMGYH0b494wDP2evvW1Xnqer6ePxTInLtVI7TMIy5y8HYoUNBVZ+jqrfEe79KRL49Xfc2jAPFxF3veL6qzgeWA9cDbwc+cXiHZBiGYRiGYRjGXMXEXY9R1V2qeivwYuCVInJ2+Sm1iCwVkX+KHr7tIvKtwsMXn0xdLiJ3icgOEfmkiNTHu4+IXCYiD0Rv4V0i8gdjjr9eRO4uHX9SbD9eRP5WRB4TkQ0i8pbSOVeLyBdE5NPxvJ+IyBlxTI+KyIMi8qxS/4Ui8gkR2SIiD4nItSKSxGOvEpFvi8jq+Fk2iMhz4rHrgN8CPhRDtT40lT8DwzB6i4gsjnbssfj3/U8icmI8Nu7fdxFGKSLnAy8HLo3Hv1Q+XrpHl3dPRC6JtuZhEXnNmPHUoq35uYg8EkPh+3r/TRiGcTiJc53743zqVhE5vnRMYxTVfdFO/YWISDyWiMgaEdka5ycXxP5pPP5NEXmdiDwOuAn49WivdpaPl+7V5d0TkWeKyM9EZFe0gTJm3K+Jc7QdIvIVEVne0y/KmNOYuJsmVPX7wGbCJKfMqth+NHAscAWgpeMvB34XOBU4A3jnPm7xQLz2QuDdwKdFZBmAiLwIuBp4BbAA+H1gWxSRXwL+GzgBeDrwNhH53dJ1nw/8NbAY+BHwFcLvzQnAe4CPlvreAmTAacAvA88CyqFYvwrcAywFbgA+ISKiqu8AvgVcEEO1LtjHZzQMY2bigE8SIhVOBkaADwFM9PetqjcDfwPcEI8/f6KbicizgYuBZwKnA2PDs95PsJdPJNijE4CrDvrTGYYx4xGRpwHvA/4YWAZsAj47ptvzgF8Bzon9ivnO64HnEGzGk4AXjHcPVb0beCPwnWivFk1iXEuBvyXM35YS5mtPLR1/AWHu94eEueC3gP834Qc2jH1g4m56eRhYMqatRTBCy1W1parfUtWyuPuQqj6oqtuB64CXjndhVf2Cqj6sql5VPwfcBzwlHn4dYeL0Aw3cr6qbCAbuaFV9j6o24/qXjwEvKV36W6r6FVXNgC8QDM/1qtoiGM0VIrJIRI4lGMa3qeqQqj4KrBtzrU2q+jFVzQlCcBlB0BqGMYtR1W2q+reqOqyqewi26n/38JZ/DHxSVX+qqkOEh1cAxCfxrwcuVNXtcTzvpdsWGYYx93g58Jeq+kNVbQCXEzxsK0p9rlfVnar6c+AbBDEHwaZ8QFU3q+oOwnKaqeK5wF2q+sU4d/oz4Bel428A3qeqd8e51nuBJ5r3zjhY0sM9gCOME4DtY9puJExM/jVGB9ysqmWj8mBpfxNwPOMgIq8ALgJWxKZ5hCdEACcRnhSNZTlwfBFWEEkIT40KHintjwBbozgr3hf3Oh6oAFvi54Dw8KA8/rYxU9Xh2G/eeJ/HMIzZg4j0Ex7mPJvg5QeYLyJJyV5MJccDd5TebyrtHw30A3eUbJEQbJthGHOX44EfFm9UdVBEthHmXhtjc1lUDdOZgxxP93ylvD8V42pfT1VVRMrXXw58QETWlNqEMO6ybTOMSWHibpoQkV8h/KF+mxCeCEB8qrwKWCUivwR8Q0R+oKr/FrucVLrMyQTv39hrLyd43J5OCBXIReTHdGK6HySEdY7lQWCDqp5+SB+uc60GsDQ+eTpQdOIuhmHMUFYBZwK/qqq/EJEnEsK4Cxs00d/3eMeHCSKt4DhCCDvAFva2jQVbCQ+efklVH5rc8A3DmAM8TBBKAIjIAHAUMBk7sAU4sfT+pH11ZHx7NcTe9qp87fb1YnRB+foPAtep6t9MYpyGMSEWltljRGSBiDyPEML4aVX9yZjjz4tJBQTYDeRxK3iziJwoIksIMdmfG+c2AwRj81i85quBs0vHPw5cLCJPlsBpURB+H9gtIm8Xkb64oPjsKEQPCFXdAvwrsCZ+Zicip4rIZEOzHgGmrOaVYRg9pSIi9WIjeOtGgJ3RVr1rTP+J/r7HO/5j4GXRLj2b7jDPzwOvEpHHR69h+36q6gkPu9aJyDEAInLCmLXEhmHMfsbaoc8DrxaRJ4pIjRDe+D1V3TiJa30eeGu0FYsIGc73xSPAiSJSLbX9GPhDEemXkAjqtaVj/wz8koj8YUzQ8ha6xd9NwOXxAX+RnO5FkxizYYyLibve8SUR2UN4IvMOYC3w6nH6nQ58DRgEvgN8WFW/WTr+GYJoWh+3vWpBqepdwJp4/iPAE4D/LB3/AmENzGeAPcA/AEtiuNTzCTHnGwhPvD9OSMpyMLwCqAJ3ATuALxLW1U2GDwAvjJmiPniQ9zcMY3q4jSDmim0R0EewId8Fvjym/0R/358AHi8ha/A/xLa3EuzTTsJamqIdVf0XwrqVrwP3x9cyb4/t3xWR3QQbe+bBfVTDMGYoY+3QbwFXEpKXbCFELE12re3HCHOt/yFEHdxGSBA3Xlj514E7gV+IyNbYtg5oEuZgtxCSRAGgqluBFxHW8W0jzPvKc7S/JySB+my0Vz8l5DAwjINCunN3GDMJEdkIvE5Vv3a4x2IYhmEYhnEkIKFU002qaklNjFmHee4MwzAMwzCMI5a4NOW5IpKKyAmEUO+/P9zjMoyDoSfiTkSeLSL3SCgkeVkv7mEYhnEwmH0yDGMmYrbpsCKEGsE7CGGZd2O1MY1ZypSHZYpIAtxLKC67GfgB8NK4LswwDOOwYfbJMIyZiNkmwzCmil547p4C3K+q61W1ScgSeV4P7mMYhnGgmH0yDGMmYrbJMIwpoRfi7gS6iz9ujm2GYRiHG7NPhmHMRMw2GYYxJfSiiLmM07ZX7KeInA+cDzAwMPDks846qwdD6WZn4y5yHLmGLdPOfq6CV0FVyFU4e9HxPR+PcWRwxx13bFXVoyfb/3d/Z0C3bR8v+3Lpmv/T+IqqPvuQB3fkMaF96rJN/fLkM0+rTMe4DOOw8MP/aU6pfTLbdNAc8NwpIXlyPwvimVIcL11NSlc90H1AOm90vD5jRyzSGXDXdUofqDTWsX26DXH3pXUf1yu/1/0cG/c6E/Tf63rj9tUJjo/fV8b0FdEx70P/7q9dO5dqf4Xa7t91HI1t3dcRKe0X7VLsdx/rtJevPfaczrnl48XnGXuv7mvFLX4Z4X1nf9ODGVu35+P9XYyL2aYOvRB3m4GTSu9PBB4e20lVbwZuBjj33HP19ttv78FQOnzg7mdQc7/CYF7nkdYCvApDeY2RvMJjo/No5QlDrSojrZRGq8LORkqeOXzmoJEgLUEyYcNbV/V0nMbcQ0Q2HUj/rdsz/uvL+39gWz9+w9JDGtSRy4T2qWybnnROTb/95cmWajSM2Uf/8Zum1D6ZbTpoDnjutECW6K+6ZyBJAuLASdyPr04gvkcckriwX7Q5ia8OdeG13e4cKgJJ6NO1n8T3ToJYEgEH6kK7OtrH1AWF0Xkdsx/P7+wX7fE6QrxO51inTfZqY8x1uu7D3tcJ+7rXOe3rxeN730e7joWttN9uVyTuiwuCC1GcaHiNbc758GMSRURJilcXXp0oifMIkDrfbksl7KfO40RxKKnLw3487sRTkXi8ve+pSKdfRXKSog0lKfVJUCqS4USpSobDk0hoK44l4jv7pVfXbu9+dUBFIAEqIiQIFXE4HIkIKQm//uyHDugPyGxTh16Iux8Ap4vIKcBDhAKSL+vBfQ6Itz7ua1z70+fR0oSRvMJIXmUkr9D0CY08ZbBZRVXw3pFlrvNHWkZh5do1e7Wvv9AEnzF1KOD3fmBrTA0HZJ8UpaX796IaxpGE2aeeMSPnToYxWzDb1GHKxZ2qZiJyAfAVgij/S1W9c6rvc6Bc+KMXkzOfPa06e7IaAMNZlcw7RloVVIWRZoXcO0RAvYu/KaXwAQf4eEEvIIo6WLlmLZJ3vO6awAMXXzTtn9GYG5ig6B0HY59y+2dhGG3MPvWGmTp3MozZgtmmDr3w3KGqtwG39eLaB8vurI+hrErTJ4zmFVp5gkfIvCPzjmaWAJBljjxLyDOHZg6KcF9PEHtSvHZ79trCLrrsV65ZixTnRFH4wKUm+IzJYU+feseB2CcFWuon7GcYRxJmn3rDTJw7GcZswmxToCfibiayrdHPaB48cyNZBR8Tp2R5Qpa7IPa8kGfhFS8dz51n71W4Rfy1Ak7JU3DN9upVutbGRoF36o1ru9rEvHzGOCjQwgTFTCCEeRiGUWD2yTCMmYjZpg5HjLjb06rTyFJyFZpZ+Ni5D567PHdkWYLPHT6XkEQlc21vnaiA03LeooCA5PHV01loWzrezuIU90XjNeNxdcrKNWtZv8oEnhFQIFd7+jQTUFWa9rMwjDZmnwzDmImYbepwxIi74VYIxWy0UhTQWPZAFVrNFK9R1BVeuzx47EQFcjqqrRBv5Uf6MVxTXYzLlJh5SSV451znl60t/jJBk+j9Q1m5dg3rL7LELEaMG7fQghlBCMs83KMwjJmD2SfDMGYiZps6HDHibnC0FvRYkQ0TyLMEcZ48d2ge3G4aPWsAEkMyJZeOV87TXfq9WIMXE6+IhooibU9eFHYaz5Go/6iUBF+qKHDq6rUWommEBw5mn2YMFuRhGB3MPhmGMRMx29ThiBF3I6OVmAUzlDtoi7g8CR47YtHIXNpJVNQpooKmCj4KNuJ52gmxVAfiFU0I9fDGhGiKjwqw1FaEaXZqokDeZ9NIA0DI96rDYRwOPEJT3cQdDeOIweyTYRgzEbNNBUeMuMtbCeJAfQzJzGP4pdO2oNO2GtOgxYoyCEV2zDG/NG2BBu31d+3QzHgdTeIFygU2Y5+2F6+4TaKs+PBq8LDxgot7+G0YM5kQCmgGaqbg7Z+FYbQx+2QYxkzEbFOHI0bc+WbSEWKeTjbM+FReWtIJsXTxWD7ml6Qoa1B6bZc7iBShnCEkU5AsiLaivwqhgo0qWmThLOvG4rhxxKJgT59mCOGfhXnuDKPA7JNhGDMRs00djghxt/xjNxBca+GlEG0Sk6YU+22vm2o4VlDsFgLNARL6aFH9IAo8dYqUf7mclmrj0fYKtouio+3zkLAgFGD5R29k0xsumdovwpgVTIWgEJFnAx8gPCr4uKpeP+Z4Dfgr4MnANuDFqrqxdPxk4C7galVdfUiDmcUoQhMTd4ZRYA88DMOYiZht6nBEiLt2mGQesl9qFGaSyRgvXEyg4qNi89IOrSwEXXnRXPt3qCiZUHjsEu1kxSwLu/Z4aCdmUaedBC2lrJoAy2++EcmEjW+yEM0jCUXID0FQiEgC/AXwTGAz8AMRuVVV7yp1ey2wQ1VPE5GXAO8HXlw6vg74l4MexBxBAW//LAyjzaHaJ8MwjF5gtqnDESHuJJe2oArCTbqLjHvayVIkCwlU2p68Un267hp2RUrNmFAlk/a+OtolECSWPOicR1s0tt/vNeD46hRNzYt3pDEFT5+eAtyvqusBROSzwHkET1zBecDVcf+LwIdERFRVReQFwHpg6FAGMRcInjuLkzaMAns6bhjGTMRsU4cjQtx11anTkq6Kgq79uyChLEFoi6/tRCjaCeN0pXMimiiIxAQp2hZ0mnZq37UfKBT75YQs0EnwUmTyFCBR8MGLt+l8E3hHBkI+sYFaKiK3l97frKo3x/0TgAdLxzYDvzrm/HYfVc1EZBdwlIiMAG8neP3MZazmuTOMbiZlnwzDMKYZs00FR4a4Iwq66KErPGvtY2PbxnrXoDsZSmzXRNs18NAo8Chdo3gt2ovDUdhpIeTKhdEzCccT7YRvxutbmOaRgQKtib1FW1X13H0cG88fPLb6y776vBtYp6qDIrYw2RZoG0Y3k7RPhmEY00qvbZOInAPcBMwDNgIvV9Xd8djlhOUuOfAWVf3KOOc/DVgNVIE7gNeqataLsR4R4q69li4KqSDICsEVk6uIdous6LkrRKEmtL1o7VDNrlIJ8dyy9w1iCYROCKaodERd0acsKos1ekVClyIBS+pDuYYUVnx4NRvfdDEr/nwNG/901RR9S8ZMQVVo6SEZqM3ASaX3JwIP76PPZhFJgYXAdoKH74UicgOwCPAiMqqqHzqUAc1WFKGlR4SZNIxJMQX2yTAMY8qZBtv0ceBiVf13EXkNcAlwpYg8HngJ8EvA8cDXROQMVc2LE0XEAbcAT1fVe0XkPcArgU/0YqBHxqylEHbtLJW6lx9DvLTDL9sZMbu8aoXqKrUVYqxcDqHw0hUCcmz4pmi4V9nL54PoK4djttf7FRQeQkLI5/Kbbgza35hzBG/RIYUW/AA4XUROAR4iGJ2XjelzK8GwfAd4IfB1VVXgt4oOInI1MHikCjuIa+5sImsYbabAPhmGYUw502CbzgT+I+5/FfgKcCUhh8FnVbUBbBCR+wm5D75TOvcooKGq95bOvxwTd4eGtr1t3Y3tyEsPImEy11WuQDvndRUth26PW7k/xSXCtYrsmyra9txJHgXeON47LbyIY7O4lPvG39/lH7+BTa+79JC+G2NmcajeoriG7gKC4UmAv1TVO+OTottV9VaCQfnraIS2EwSgMQ7eJrKG0ca82YZhzEQmYZv2l6tgMvwU+H3gH4EX0YmQOgH4bqnf5thWZitQEZFzVfV2wkP1k+gRc95Cn/Jna9qiS8oCTKVb6BV9fHxTCs1se/GgXdsO6KyLK5dTKIu29vklL5wvvIRjEqeMh+7jQHmu6SWsxcuFjf/X1uLNFfJ9/ewniareBtw2pu2q0v4owTjt7xpXH9Ig5gDePHeGsReHap8MwzB6wQS2aX+5CgAQka8Bx41z6B3Aa4APishVhOinZnHaOP274gNjJvKXAOtineF/BXqy3g6OAHGXjghZXyfzpXg6HrvxUkwUnjoFLXn1cKXQzYJcOs61nFi6QMPauKp2wjWL6zpKwo8oOKXtqSs8e11r/8ZSEonlMZuwmzvYk/GZhWXLNIwOZp8Mw5iJTIVtUtVnTNDlWQAicgbwe7FtMnkOUNXvEJe+iMizgDMOabD7YcJvQUT+Enge8Kiqnh3blgCfA1YQMsb8sarukJBe7wPAc4Fh4FWq+sPeDH1y5HXFZdIubdAVilmKwIRSjpXxLjTGw1YULB+XUrhlKJEQ28vJWGIdPS159/ZK0DJuDbwgHsf2XX7TjWx6o5VKmAvYmpbJ02v7pGqeO8MoY/Zpcsz2uZNhzDZ6bZtE5BhVfTQmR3knIXMmBC/eZ0RkLSGhyunA9/dzfo1Qcuq6Xo11MhL3U8CHgL8qtV0G/JuqXi8il8X3bweeQ/hQpxOy7n2EvetrTSu+pkgLJKcUktndp+woE4Jw82mpf0loSTl5SvTAFQXLyxcsZ8gMrxoW9RWaLCZUaRdLL5+uAnlp7d14oy3dQ3wQeys+sto8eHOA8PTJBMUk+RQ9tE/mpTCMbsw+TZpPMYvnToYx25gG2/RSEXlz3P874JMAMafB54G7CKGWby4yZYrIbcDrVPVh4BIReR4hdu8jqvr1Xg10wlmLqv6HiKwY03we8Ntx/xbgmwQDdR7wVzHr3ndFZJGILFPVLVM14AMlGXJIDtmABn2lWmRO6VBKmhJE097t5ULnxbG216+dhZPudXRjhVlaWmcXC5trWUCWxlQWbmXxJ3lH3LXLOpQSrBizH0s1Pnl6bZ/MS2EY3Zh9mhyzfe5kGLONXtsmVf0AwcM+3rHrGMcTp6rPLe1fQiif0HMO9pH0sYXRUdUtInJMbD8BeLDUr8gYs5eBEpHzgfMBTj755IMcxv45+5J1VGtB/KgDTTWKtBgCOU62SxVCTbuxMZpjBVhxerk9rs3TogRC0ddp9zWK/dICvqIY+liPn8axFiUUgK7wzq6+5aLrxqxFgdzWeR0Kh2SfyrZp8fE1m8gaRole2qd9hS2O6bOc8NQ8ASrAn6vqTcwOpnTuVKe/t6M1jFmEzZ06THW80YQZY9qNIf3ozQDnnntuT1TJ8DKPX9qi754avqb4Po9kgoyW1quNHWUh7LSTLKW9bq7w6kVPmrRrz9HtuSt+t1y8gB9T885p2DIXwjoz6RZ+Y7x+inbXwVMJ6/qks1avEHYrPrI6JFh5k4VnzlYs7KlnTMo+lW3TyWcvUEuoYhgdemyf9hW2WGYL8Buq2hCRecBPReTWGPY0WzmoudMCWWJPdA0jYnOnDgcr7h4pQgZEZBnwaGyfVMaY6cL3KTqaMLIsR1NFWoLkgu9T3LAEh1spk2W7lh2EzJdlL11RdLycLjOeq6l29ouyB20PnI45J7YXgi+GXYqX4JFL6PYoRoGoqghhLV67UHrMtjkep3xwDRvesurgvjjjsGOhgIfElNkn+2dhGHvTQ/u0r7DFNqraLL2tMbsWJcyKuZNhzFZs7hQ4WHF3K/BK4Pr4+o+l9gtE5LOExcC7DmfMuLQEsiCiJJOwefDVULKgKCgOBAFVeO0KD13hrWt73OLxYm1eomgliMa96tXFNXVhINr972esQCR64GJ7e52ddPdVUUQ6IlBU2uUWXOaCyGx79w79+zMODyYoDpkps0+q2M/CMEr02D7tK2yxCxE5Cfhn4DTgklnktZsVcyfDmI3Y3KnDZEoh/D/Ck7SlIrIZeBfBMH1eRF4L/JxOMeTbCKl87yek8311D8Y8eTQmSXEd8aUO3KjD93tcM/4SjElmUq5D1ymXIO1adu1rVYOrT/tyaLpOjbr2tUqir+xxk7JwUyjW28UTVTQWS4/hnz5e13dCRMsC0DVdV9KVjRdczCkfXAPAadev5f7LLpr679boGYrVVpssvbZPilgMv2GUmIR9Wioit5fe3xxDCQEmKhI8uTGoPgj8LxE5HvgHEfmiqj4y2fOng1k9dzKMWYjNnTpMJlvmS/dx6Onj9FXgzeP0PTwUiSVzujxrvhbFUxISrIQsmaXwyTGJSorMmEGkdbJXSsOhNd9eb6daymoZr1UItEKchTV9BI8i4TzSGA8aC6wXa+lIFBl1pbV1rl1DT4tzimEnQN7JplmEZJqwm31YNrrJ02v7pJjnzjDKTMI+bVXVc/d9/r6LBIvIvsIW93Wth0XkTkJh4C9OMPRpZVbPnQxjFmJzpw5zu4BTeZmcdEodSCZBQ3mJZQ72IexKmTJ/ftXV+D2DALj58zjpuqs7CVNaAlUfXstZOOl43bQQj0WdPJVwfiHSsnhCu16eIA0HScmTV/Od69ZzdCj8+LSeIy2HVoInb+W6NfiqsvHNllRlNmKCYmbhx811YBhHJj22T/sKW2wjIicC21R1REQWA08F1vZqQIZhzA5s7tRhbou7sRRzNA+0yyKUj8ewyFJNO3UhkUkh7ICufZLgYSOTtkAUjaULijBN6A77dASlmYSEKhDFWtMhVY+OJJ3Qz1J5A8kkhoKCNhOoaNvbpzWPNB2aKD4NIvOUD6xBPKy/0BKrzC4sFHCmYGGZhjGWnv5NjBu2KCLnAm9U1dcBjwPWiLTjaVar6k96NSDDMGYL9v+6YE6LuxCOGTxkWvLEidIJiywSn5Tq2oknJFeJHrZ9JKREE8XVM/xo2hZpAFrxSMu1r9dO2lLggWpH2HXq4YE2XfACNjrrAcVL8MolsYSCD6UQxGnoX/IwSlEmAfCV8GqZM2cX9vRp5mALtA2jm17aJ1Xdxvhhi7cDr4v7XwX+V08GYBjGrMXmTh3mtLgL6+m0lCBF8FVF8tIaOmJyEy+dZCkQBBgdzefmz+sKy4QguvxQJQqu6LFLfViXF08UlXDvGBLqhhL8ggxUkIpH4xo5bYW1dcELGIShpgq1HG0kob0lITlM1UMmqAZvoeQSdGhRj49YO49iTBq8eJmwfpWtwZvpKILf1xMFY1pRhczbk0DDKDD7ZBjGTMRsU4c5Le5cFDh5n8eNuhAJ2YieuHaZAmIIpHZqzxXJLEui7ORrrwbonCsdIYWX4F0rsmbm0ilLUKzxi4XL/fwMWtE7V3juopcvmdfCNxOoeFQdOCWp5fjdFYjX03q4B6KhDUVdDP9sSbtenibaTq5SfE5f96xcs5akCT4Fl8F9l5vYm2lY+v2ZhdXNMYwOZp8Mw5iJmG3qMKfFnU+VpCkl0VZ4swApkqrQXXi8nIQF9lseVUuhnNqfhzDJigefQC2HVsym2XDtDJkykKO5kNYzstEUckGaDuYHYacKEjNlumqO31brSqRCVkre0nBB4FV88Pa5WE7BRc9dzO7Z9kjG76G1QHGN8ITj9OvWct87TODNNOzp08xAETJv/ywMo4zZJ8MwZiJmmwJzWtwljfBDTkdC2OJ4Se86mTJLqs51RBvxvI7HjnZGTFnYCt2GU8QpknhcBdxAi6yVhNIITUeyoIXPBW05NHNUFjZo7aqFi1Y8zMvQkbRTB6+WQxHyWfcdL2MrrK8Tp2jucE2H72uRVD1+TxrW5REzdBZex6RUigHQNHoSk87awjOvXsc9V184VV+7cYjYOq+Zg2LZMg2jjNknwzBmImabOsxpcQd0hNmY91KEYkLHqxfDLdv9C0FUiC5fro8A2kiQig9hlF5wVcU3E3KfUB1o4lNHZUFGs5mSVjyt0QRSJWukSD0PY2kkaLvmnSIVT1rJyRSqfS28F1q7a+CFdGGTfHstrLWrKIwCuZB4foriAAAgAElEQVTvqUCfjwlgwjg1UUg6de+0/FkEpFWEigp5v3Lmu9cBcM+7TOQdbsxbNHNQFVr2szCMNmafDMOYiZht6jD3xd04yNgMlmXNFmviFfvtJCVKJ7yx8OA5RaMXr2/hKHnuyBsVJHMkC4Lgy3NHvd5iZKiK68vwIyk4QRKPJEredCDRE5goIor3Dh1JcQNNGoO1sPZuICPfVgulELyQDLTQeobbVgsJWkZjiKZKW4xKS4KXLi1lAxXAg695XMy0KRntRDNnvnsd6uDeK03kHU7MWzRzyCy1smF0YfbJMIyZiNmmwJwVd6e9f20QZ9BdY67sySvtS1HjrvBwOYKYc7GIeNb5hZHCM9ZwUPO4ek6rmeLzsL6vsnSEJy57iNt/fjLOKd4LfiSlsqARQi0lePhoOmQgg12VEFK5Jw33X9iksrDB6GANV8nxrQSfhXIIbiQhHRSyRh1f91RGBT8QwjClFgubVz1uTxoFnOArijRDH614pBFCOrWiuKaQ18NnT0YEF3XgWe9ax8/ebQLvcKCKeYtmCArkFsNvGG3MPhmGMRMx29Rhzoq7QtiV8qgECkdcex0dXWvS2p67UnKV9hq2IsQx1fZ6Nkk9zil5y7Fw0TBDqWfZ4t388KETQZQl84d4ZNtCUGjtrCN9GdpMSLel5Mc30MEU+nLIhfqyYUa29qOZww9WqCweJWukpPUMnwvVR1IaJzbxzQq+6pF6TmtRFHWZCyUT6h4ZTNp18SQL2TMlF3xfHj5QEhKvaKLkfZ5k2OFr4YvRVEPeFQ9nXbUu1AoE7r7OhN50Yel8Zw4W5mEY3Zh9MgxjJmK2qcOcFXdFApV2OQOIXjg6Rc3H+x0o1q2JtOvXEa+FC+UF1INWFemLGTJV6JvXoJkliFMe3r6Ap67YwIbdR/GLnfPJGwnp/BYu8bRGUyrzmvj+jOoDddRB8zhPZWGDxfOG8XfNp3GcIi3B5wn1gSaNkQo+czSOa0HLkS3Kgodwc5XmsVlY+9dwoUzCaPDIVXYmII7W0hbJjkpbrLrBBK0qmna8eeF7UvI6uFbIIpo0Q7IVTULY5uPeGcI1f/YeE3m9RrFQwJmEhXkYRgezT4ZhzETMNnWYs+IOaK8vAzrCTrTbYxdpt5UzYyqdRCsl7554gYZAViE5ahSA/lqLrY/NZ8HiYfprTX66dRm79vQxf94IjZ118pYj21NBRhNa9RxpOTSBvK4hrHO4ypbGIvSkJmQOTZV8T4WfvfQKln/0RhDYdP4lnPKZ9+IHK0guNI/NqDyWkteDWMMBaVifl7Wkk0wl1VDQfdShNQ3r9kTBJ0gzJFZxzfDBfaqIp712TxOgCq4RSiuc9a51oCbyeo03AzUjUBUrYm4YYzD7ZBjGTMRsU2DOirukEdaa+UqnLWg0QbyWRFzppPIaPB+zSYrudVxiSKamwWu3eOEQgyM1lizdw9L+YRbXh3lgx1IG+hss7R9mZ/8AtXqL/OfzcS2heXoDBPLFkDcSGKxALUdHE9KdKdnSFm5+hnuozimfeS8Dxza56wVX84RV60iX1OGMIXT9ALUdCUNnj5I+XMN78EMJ9WVDtB6Yj87P0aridqfgwveRLQolFsilI2RjxJlrCb6mIblKJXjrxEMyKu3PDqH4ufgg8iSzcM1eoCr29GmGoGBhmYZRwuyTYRgzEbNNHebkt3DGNevI64qvaqdmHWMyZBa6TbvPbdclb9e6i6GLSaeD789DYfFEyfZU2DXYx8i2PnbuGuDee4/njk0nI6IM1Jrcu/E4Tlv2GGmak81X0l/eyYaXX4EfrOC3VSEXaktGmLdoBHIhW5whowmVDXWyBTl9/U2ajZQVn34fu8/MuO8dF5E91E82z3Pn9Rey7JidJKNCPs9T3eEYHa6SLczBgRtx6JIWKOQDHhlJQgkGD9Jw4VVBa9pOquJrStKQIOwaQfAVnr+QObQ78cxZV66b2h+eEQWF2+82ESLybBG5R0TuF5HLxjleE5HPxePfE5EVsf2ZInKHiPwkvj5tqj/fbMMjttk2Z7cDZSL7ZBiGcTgw29RhTnru8nohRmgX/S7CKtuirRSCKb60Xy6J4OO6uxjCqYm216hJ05EsaZC3XCiHkCi6s0rtmGHyLGF+rcE3nraGZ33zbdz78+PQhuMJv7yRex85muU33wjQLjre/EU/jaqSLGqS76oysCnhzuv39oj95lcvBeDkJ2zh35++GoCHNy+hL4dkfou+ZbtpPbAIXZjh9qRoVWEwxc/LcYMJCPg8lEdIhx2tJRmMJLhhB6K4TMgGQqKVyh6Hr4TvrR3emdAO15RMgofPB4En3rx4U8mhLAoWkQT4C+CZwGbgByJyq6reVer2WmCHqp4mIi8B3g+8GNgKPF9VHxaRs4GvACcc9GBmORaWaRh7Y0kLDMOYiZhtCkw4axGRk0TkGyJyt4jcKSJvje1LROSrInJffF0c20VEPhg9Av8jIk/q9YcYS2V3KZTQhaQq2l5zR7ew05h0BZC8k3AlNChaUTQqPskFib846VGj5I/VEafMnzdCZaBFsqRBY6hK3nI8NjjA73/rAhp5ysZXXIbUPBt3LKY5UqG6OKzTqz2WsOm1l6J9nmRBk3x3lU1vuIQ7r7+QFR9a3f48p65eC8BDdx3L8o/e2BZ2ANJIGDkhIx9OGf3hEiqDEurl5UGMas2T7E6o7nBUtzuSvgxUSEaE/o0VatsdLo+JVBwkQ+HLaM3X9veSNKIAjt+JZBKyaApoSjtk07x4U4MSQgv2t03AU4D7VXW9qjaBzwLnjelzHnBL3P8i8HQREVX9kao+HNvvBOoiUpuijzalTIdtUsI/C9tsm6vbgTKRfTICs3HuZBizGbNNHSbzaTNglao+Dvg14M0i8njgMuDfVPV04N/ie4DnAKfH7XzgI1M+6gloze+EY+oYQdfeKHvpFHUxKQnRU+VCIhLJoqDTUCNOE4U0Fhpf0MKPpLTyhEULhpGNfeGeLUe9knHXlmM5YWAXv/Gvb0dST5472FUh39zPpvMvIa93xuwfq5MMdn4cGy+4uL3/wMUXAbDhravY9IZLWP7RGznnS1eGfv839Ks8lpLXlOaKBq4a6hdIFtfXeSHrD8lR/LZaOxxVkyDaXEPw1eCxc61wnmuCr4bvqrlAyWtaqh9B+H6SsEYPolfPwZlXm8A7ZHRiQQEsFZHbS9v5pSucADxYer+Zvb1v7T6qmgG7gKPG9Pkj4Eeq2pjKjzeFTIttUhXbbJuz2wEzgX0y2sy6uZNhzGrMNrWZMCxTVbcAW+L+HhG5mzAxPA/47djtFuCbwNtj+1+pqgLfFZFFIrIsXmdaeODSi1i5bk1YU1auXwfdNe2IYZoqobi3CmQhXFKTUI6gCOUU4to7FHEhk2Uy5Jh/2k6Gf7qYhsL9V1zEKZ95L5oIJy/Ywe3PeS8AZ3zxPWx42RUhOcpJe8h+tIjTrl/LwBbhtM9fA64aBOPCoJTOuGYd91657xDHTW+4pOt9uqBJ1qiTHjtMsmmABY/fxq6H+trqtbo71LobWdEi3ZZS3ybktVDuIK+GbslI+KwuF/LEoyq4FvgE0lHBJ0EAS7vGBKDBg6cpISupgrjgwfvZNRaiebAUceMTsFVVz93HsX0V+Zh0HxH5JUKo5rMmGsjhYjpsU/Ek0DCMwCTt0xHPbJw7GcZspte2SUTOAW4C5gEbgZer6m4ROYoQAfUrwKdU9YJ9nL8E+BywIp7/x6q6oxdjPaA1dzHpwi8D3wOOLYyOqm4RkWNit315DboMVPQ0nA9w8sknH8TQ98/6C1cFgVcKuywnSSm8dp2smEAOpFHYNQWthLp2IfzQt4uYS92jqSDLmow2K3DaEPe96EpWfHg11WMy7n3hVV1jaTVTln/q/STbazQbdVorG5xy4mPsGO4jf2Ax1WHhvssv4pL/fhEAtXMO7GedDVVIG4K7Zx7NRZ6d9y3BL8ogF5IhR14LXrvqIyl9jwrpSKjfN7xMSBrgBmH0aEiGJXjrCtEWv7e83vHQBU9mOCZFQpYomMWHZYpaU864di0qcN87LjqwH5wxFYU4NwMnld6fCDy8jz6bRSQFFgLbAUTkRODvgVeo6gOHMpDpole2qXbMfHKbyBpGGysUfOD0yj7V6e/puA1jNjENtunjwMWq+u8i8hrgEuBKYDS+nh23fVF47a+Pie4uIzzYmXImPWsRkXnA3wJvU9Xd++s6TttYrwGqerOqnquq5x599NGTHcYBsf7CVay/cFW4XxR24qUTlthOqCLBS1d8gKZ02oTwmse6cZXg1ZLhhDxztH4+gIvX3vimi/cSdgAbXnYFtARN4b7LL4KhlG88bQ07tixgw1tXhTbgxnO+wBNuvYrROxcd0Ofc9NpL0UQZeNJWWNiietIQye4EPFR3OtRBa2FOa4GnuQCyuuByWPadBvMf9NS3K/0PC7VdkI5AbXvhmYs3KAlkX1HyqpI0wTWjyCvW37WT2Ejo7+CMa9ce0GcxArm6/W4T8APgdBE5RUSqwEuAW8f0uRV4Zdx/IfB1VVURWQT8M3C5qv7nFH6kntFL25Qu7D/sa6Jss62X28FwCLbpiKOX9qnCjFwObRiHjR7bpjOB/4j7XyUsXUFVh1T12wSRtz/KuQ5uAV4wFYMaj0l9WhGpEIzT36jq38XmR0RkWTy+DHg0tk/GazBtnPLBNZ0QQl+ynTGZSnvdnUZBk8ci5cU5rRCOKS0Xin3XQq04d1STjX9yOZIL9/xREHTL//L9e93/6d8Iwm3T6y9lw1tWccpn3tteJ7fp9SH75bn/ckW7/09+/z1tsXcg5PNzfvR716GjCc3NA+gxDagoo0fn+AokQwkDDyb0b1FcpjTnC1lfQtJU+rZl+JSYGRN8LSRYaf+rkU7oZtIQqrulvfROk1JdvDwmpUmjx9OHr/X095nAOxBUIfduv9v+z9cMuICQ6fJu4POqeqeIvEdEfj92+wRwlIjcD1xEZ93HBcBpwJUi8uO4HcMMZTpsU65im21zdjtQJrJPRofZPHcyjNnGJGzT/nIVTIafAsUc6kV0/71Ohi6vPdCzudVksmUKYSJ4t6qWZ+nlJ/+vBP6x1P6KmPnp14BdhzNmfMNbVrHhbatwWdgvRAg+bhRJV0p17QrBV5RQyEKpA+o5fiTlN854gLzpOP29a9vJTgA2vWZv7+q//U63sNnwsis4+5LupCPF2rwDZeWazrUr21JOu2Et8+9Lqe50+KFK8N4lwPGj9D0mZH0hFHPXGUo2AFmfkA7lZHVHOqIMPJKHrJgesv4QqiqE964laCV8d1l/SLbiK9Fzp0EQ+kopCU0ev9KYpbTI+GlMhkN/2q6qt6nqGap6qqpeF9uuUtVb4/6oqr5IVU9T1aeo6vrYfq2qDqjqE0vbo/u71+FiOmyTTuJnYZtts3k7iL+8Kb7e3GS2z50MY/YxoW3aWni943bzXlcQ+ZqI/HSc7TzgNYTESHcA84Hm9H6+yTOZNXdPBf4P8BMR+XFsuwK4Hvi8iLwW+DlBxQLcBjwXuB8YBl49pSM+SIrwzPUXreKUP1vTXiNWrBdrr8MrQjarQdz4mo8lFRQZTtGq57/uW4nbUeG+K/btYVt+y/VseuVetaMB+OmNU5NsZP2qzv3FQ2uBZ7giJMOhlt36C1ex/BM3kO+p4FPI+5TatuB9ay5Sdp+SMG+zkDSUhRtbDC6rkDRg+PioerMgdJMMfKq4prRrCLpM2sIOX4onyQRNgsjTpFSCQuGUP1vDhretmpLPPpdRsCfgk6P3tkntZ2EYZcw+TZo5MXcyjNnCVNgmVX3GBF2eBSAiZwC/d4CXf6RIkjTGaz/lTCZb5rcZPxYc4Onj9FfgzYc4rp6y4W2rWPHna9jw1lWcunotkkG2wAcxA8Gf6UPWTBSoeEiU2pJR7vmjq3j8FetY+rT9R0tseuVl+xV4U8kZ16zjviu7hebyj4ZC6elAi2x3lbuvvZAzrlnHnTdcGOrROcjqMHiCY9EDObtPqjB0QhC4Po3ezbieLhmVkGiFkF3TZeATJWkIvhKrTuQxoUpCLDchoWwE4PKYdXNEWLl2Db4C1V3Cve+0ZCvjUiSpMfbLdNim8M/CvBGG0cbs06SYi3Mnw5jR9Ng2icgxqvqoiDjgnYTMmQdC4bW/nm6v/ZRzQNky5xIb/zR4kIqwypVr13SKmftQ685LyA6Zbk/IV44wuqfGmVev45737tvztvzmG9l0fihV0Cthd86XrmT3nr6QqAX2Kpuw4kOrYX7OGdesIzshxY06znrXOvIB5ZwvXcnPrrmGsy9Z165zt3NlgnhoHJOHjKCi4ARyqAwFL59rEdbkxeQqriXktU5pBI31A5MWZEnIsEk4FMSij3X18lBkvRCLB8oTLgohrT9ZO3dLLShhUbAxMzioWmCGMUcx+2QYxkxkGmzTS0WkeADzd8AniwMishFYAFRF5AXAs1T1LhH5OHCTqt7Ovr32U84RK+7Gsv6iIPbOuHYt9ceE4WVh7Zg6aC3yOBXSrRXuuXr/oqIQdlPJKZ95LzxWw/d5pD9j459cA8BvfvVSNv/s2LZQLSgXQF/x4dVoRRk5qYUbTBi8ZzE8P4SGnn7dWiqDwvAJOTo/Q/akMctMOFcTpTUvvFa3OkaO1lD0PK5brAwJPg3r76QVzmnNU5JRgZEYwllk0YwF4guvYF5XVq5b0w6XnSxzWdR1sLUrMwVVsRA0w+jC7JNhGDOR3tomVf0A8IF9HFuxj/bXlfa3MY7XvheYuBvDve+8iGc89Vo2P22A477X4Jv/GrxvT/ny5exI82kfz69+5TJq9T70xKydlfOcL12Jc57tDx2N8/s/v7rdkfUrx526jcd+cGy77tzyT72ftFLBV0H7fFhPGLNcFuGUiOCGoP/njj0rfTjmQohlkkNrAFwrCLu8X6nuFHxcp6cxDLO5QDvhmZngfEfo5X0TDP4Ixlso4IzBItAMo5te2acDKfIrIgsI2YD/fl9Fgw3DOLKwuVPAxN04fO0/39n1/uyL1zFy7NE8cMn0rBF73N9fDcDdf3A13/vd63n8P1zNiiXbOPWz1+Fz4ZilLQZHayQLmuS1hFNXryUf8O0SC0/8p3fy4+ddGy521iDrY+295RtuZPnHbojFzSu4DBpLfEyCEqawGjOEas1TfSQlm6fsPl1DiQOBvseEyu4QVulTyGsxW2Ym5NUg3BrzIWlAqx7LI3hBU8VXQ1IWJWbVVGH5TTeSDCVdyWGOdIp0vsbMwNvPwjDa9Ng+HUiR32uAf+/VQAzDmF3Y3KmDibtJ8NPV0xsKePcfXA0Ekfbsk+6mVjmLe767gvz4BscevYvtuwdQL+SNhL75DZrbK/zmk+7mCbdehVfhzvOubV/r3hdexSkfWINf1GLTG4IXcuXaNfh+z/pVF3P6F64h29KPr3uoeNK+DP+LOr7f0zwuQ0YdeCEdESp7hNp2xWWQe6W1VMhrsZh5oiChlqBraafoeV2RLAjBpBFEZKiFR1zDJ+HeRheWsGBmoNiTQMMYSw/t03nAb8f9W4BvMo64E5EnA8cCXwbO7dloDMOYVdjcKWDibpKc/oVraO6o92RN3b4YHKrz2f/+FTQXZL5Htlf5xfBRoSzDqIO6Z2R7H1JRvrtpBc4p8/sb/M7XV/GNp60BQvjlpre+nZVr13Dq6rVh/dxOx88uCoJVvaCLmzinSKJkwynMy5FRF+vaCdUdDl9V8hqMHC3B47dYqeyBZDSMNR0W0hEYPUpxLSGb50OphCysuUsagq9qqC3oQqIaBHzVg4PlN93IpjdO33c7k1HEvEUzBRVLqGIYJXpsn7qK/IrIXkV+Y6a6NYQyA9OyfsUwjJmPzZ06mLibJM2h6rQJuxWffh+VWkbWSsNTiJYLhcQbArhYhF1wg45snicZdejPB2jMy2lIH9sygafFi0WnWHryEM1f9Adh956OJzJrJmgzIRloIc6TtxwykIXwzJZDEyWvK9l8j7SE2naHG4WBh4X5D2bsOSEN5RAyaC6A2g6huVBJh8KE2FeC104TRV0IxdRUQ1hm04GAG3WhnqDRxh4+zRzMc2cY3Uxgn5aKyO2l9zeXiwWLyNeA48Y57x2TvP2bgNtU9cFQJ9wwDCNgc6eAibtJsHLdGtyx05NM5cy/fQ9HLWmw9ZEFQdQ5JRlyuKbgWjEpSSPUo3NNSEYTJIdsQKk/mpLXlGREWPGh1Wh/TtqfsfyW65GkQvW4YRrzq5zxxffQ2FVn02svJal4chWyoQp9i0doCehQSmVXgk9Crb/W0gxagssceU3RRKjuhKzumP9wxsiShOb84IlTB8lIqH+naciu6ese1xBcw4XagS6uw6t50Oi9E1jxkdWQd2f7PCLR4FE1ZgYW5mEYJSa2T1tVdZ+hkvsrEiwikyny++vAb4nIm4B5hNTjg6ra+6KyhmHMXGzu1MbE3SRIRoT7Yk25XlHUx7vnj65ixaffR9KX4Udr1LYEwVbgGoLkkI6EZCXV3UpeF/q2Bs8ZKqSjoDsSGlWlr7/BqKug3uGc4iqe7MEB5KhWvOc7OO1z19IarfGzP3xXSLiCkB3fQEcTSBVaQv2R8KvSWJqTjDiSEWHXSofExXUqkPUpAzugNRCEXp6G8VZ2u7b3TrKQSVNTbXsVIYR/kkkoo2BYKOAMQRXUwjwMo4se2qcJi/yq6suLfRF5FXCuCTvDMMDmTgUm7iZB84Rmz+9RDvms9zcZ2VOD/gxfDevdkoYgmSBRELkWpMPBC5ZuV3wK6TBUBkO2S8kFTRMGWwuoHTdMsyU0RipI9Kytf0Xnf+H9Lw7ZQVd8eDWkAlWPDqaQaqh9B/g0ZMhMhh2VPcLoUiVpBlHnsnCdyqDQmg9ZPQg4NHgYJQethzILPlXwQaAioKLt9XeuFTx/K9esxWVw/9uPzAyalsRjZmGOO8Po0GP7NG6RXxE5F3hjuWaUYRhGGZs7dTBxNwnm3VXr+T1e9t3Xc//OpXz/2e8DQATcY1XyagjFTEYFyWiHY6YjQVylo56kFQQdAmgIl2xmQey15iVkj83Hn9IgqeX4zLH+ivFF08Y3dYdDnvmedYwek4ED1wp18FoLc/J5QrrHkdeVBfcLo0cJ1V0wclwsVJ5AEvVw1h+EXqiR1ymQHgrEK+IEyaNoFQANIZsCp964dtrKT8woLLRgRmE/C8Mo0UP7tK8iv6p6O7CXsFPVTwGfmui6ItIPrAJOVtXXi8jpwJmq+k+HOmbDMGYIc2zudCh2y+KNJsGd7+99KYTP/NrH2sJu9JEBkmpOMhpEkmsU6+0Ahb6tSmVYqe3OqQx6klGlujujtqNFZTCnviOjOqikwyFsMx2GvvuDQN34J5fzxH96535GAudcsI4zrl3LPVddCFXPpvMvobnYk9cVaTmkKeT1ICj3nAKt+crgybH8gScINRXywoNXzpUSPXkUYk7DpsU+tDNq+oqy4i9Wc9oNa6fse5416ASbMX1M9LOwzbbZvB0MU3293vNJoEFYrwewGbh2390Nw5iVzD7btD8O2m6Z524SLL/5RnDKptddyll/925+9ofv6un9KktHaO6owzEZbiShujOEMboMkhEFhXRUqe7KcI0cUQWvSO5xicPXUvoeU4aWVcjjGjxRaPyiDtApcL4P/vtDHTE7sGSE5bdcT31njXvedSGPv2IdjSWKJrHkQb9Sf9RR3xYSrTQXQGuhkveV/pJKwk3Kf2At2es4GkskVDSuyxN8ckhf5yxE5tTTp1mN2s/CMLqZlX8Tp6rqi0XkpQCqOiKWatMw5hiz0jbtj4O2WybuJkF5PVyvhd2Tb3sH8/phSJTGjjrSEqQFKCSjSm2Xko4qtR0t0j0N3GhGyPoQBJ5WUlwjw/dVqAylJM2wBm74OEGXtA54PHe94OpQ9DxVTr1xLflJHnWKazhQ6PuFo+9RxVcg6wvCDoJGg8KLV+zH0EwAKQm9wrOnITxT0+AFlFYRrqmcesNaHri0O0TzjGvXkowId183vUXme44yZxYFi0gCfGV/GfJmPHPkZ2EYU8LstE9NEekj/jcSkVMJT8QNw5grzE7btD8O2m6ZuJsEK/58DRv/dNW03OuO517Hio+sZvHJO2gMVcmPatHyFdgtJA0hq0NtZ04y3EJaOTI4HE70Ck4gz8E5klZOrS9l9KhKyF6pwNDB/bjXX9T57GddtS7UsRuBvArNhUo2L9Szk0zBE4RZIepcEHFKnCOXQjHbxMQrAJqE5DBAu9i5JgoCp12/lvsvu4jTr1uLpiAJc0/YFcwRA6WquYgMi8hCVd11uMdzMMyxJ4GGcejMPvv0LuDLwEki8jfAU4FXHdYRGYYx9cw+27Q/DtpumbibBBv/dBUr/vp9bPw/l0/L/apHD5M4pTbQpDFURTxUhmHelhzXUrIBR20buMd2AqAjI1CrQZ4jaQr1Gr6/RlZPSEc8KgkDWzx3vv/Q6sedunot2fE5okKzHoRX0ghhkyqE36a2Fy6IOmkS1t05OipPOn2IzRQeu7STNRMJYg8N7zVRTl29FklDJs3T37e2LfjmHLMzPnxfjAI/EZGvAkNFo6q+5fANaZLE303DMErMMvukql8VkR8Cv0b4D/RWVd16mIdlGMZUM8ts0/44FLs1YUIVEamLyPdF5L9F5E4ReXdsP0VEvici94nI50SkGttr8f398fiKg/5kM4jpEnYAxy3aw7Yd88izBFoO8UHwjCxx+FTo39IIgkgVv3sP1Gr4nbug2YJ6DUYbaCIkDU/tsVHSEeX7f3XonscHLr4I13SoU5IRF0sZBPFWeOo01SjkwrG8FsWd0F5Pp3ETL4iXtodPHaEOngsZQouyD+ridYtzXfDiSSZzV9h52f82u/hn4ErgP4A7StshMW22ScU22+budqBMZJ9mIHGdynOAJ8dMc/0i8pQe39PmToYxncxC27Q/DsVuTbaSyf4AACAASURBVCZbZgN4mqqeAzwReLaI/BrwfmCdqp4O7ABeG/u/FtihqqcB62K/Wc0pf/NeAN50x5+w/JO9/Ti/+dVLGW5V6OtvkmcOGXWkg8LAL3Kqg8rg8Qm7V9TJ+1PyE5YilRRGRnH9/WizCbmHSgW3Z5TKjlE0dUFcHSKn3riWFX++BuisnXMNCUIsh2xBHpKsRFEWPG20BVrhycN3h2CqxL4unAuxbIJ0zqPoI4ALyVw0gXuvnKMhmYD6/W+zCVW9Bfg88F1VvaXYpuDS02ObvG22zeHtIJiFtunDhIxzL43v9wB/0eN7HvFzJ8OYbmahbdofB223JhR3GhiMbytxU+BpwBdj+y3AC+L+efE98fjTZ3tWqiQNvxUffvKn2fTqt/f0Xt9+5g1s3bCExHlc6tGBnL5HlfrWFq6lVHcrrQEhGWqRz6vC0iXIwgX4wSFUFR1toP11/MJ+fD1FnbDk7hF++Y2HVk7ggUsuYuOfrgrCKvX4CqQjYT1c9v/be/Mwyer63v/1Oaeql9kZNkeW7hlmQFCj5k40PsnvdzUmxC3ivVcNxlxRJDxE3JhhGUAQkWWYpUeIRkBA0WjQGBeu8apIxDx5EhVcYpDFGaBHJqBss/Z0d3Wd87l/fL9nqerq7uqlupb5vJ7nPFV1zvec+p6pPp857/PZFsaEB7yCjF2rg2AMwlHScErA97nL+cyT0E3f0DwuZHGaEjuRh7rvkBjXHqGQHbOjWyTM9dP2JiIifwb8HBc7joi8VETunO1x58U2JcmittjSqcuMLr45Pl7jeYWqnosLEUdVdwNdjfxCu3cyjCbQfrZpMmZst+rqcycioYj8HHgKuAt4BNijqmU/ZBdwjH9/DPC4n0gZ2AscXt95tCa6cwGr7rh63r5v8L3nEwYxCxeMQilg9DBhz5puyr0BhVEngPatXkQcBujCHgiEYOVxBN3dEEdodwEplRk6toehY3o4cGwPP7txbsIXB889n51/dSEUlNKS2HnqykLcHTuPWwipt60AiKZizIVXOoEmPo8uLrpwS2LfzDzyoZoCGrh2CE7kue9Pmp0/um69E4WdiLrznGxpM64AXg7sAVDVnwMr5+LA82Kbmt2HzBZbGrlMF21L2zTmK/cqgIgcyYz9lvVzqN87Gca80p62aTJmbLfqKqiiqhHwUhFZBnwNOLnWMP9aSx6P+y9ERM4GzgY4/vjj65lG0yiu3k+xEPGHd13InuEeVizez12v3tbQ7/z5G6/ild+9iOjogIMjiyktFbp3C4Vh0FAYWQZhqcDoYUsIS4vp2jNGuND1sdNiwMhRvSx+9ADfvfeKhsxPxoTu5wJGl8doQQmH3HMCDSHqdeGWiRBzfevwbjogFmJfOAV1oZ3ibzS04PrnqX/skHjwNFQ0uRkRWPnxrUhPQ06tBWjbp0wTUVbVvVUPoWdyWzmORtumcPky97DBMAxPW9qnG3D24SgRuRp4C/DhRn9po+1TDwvmaKaG0Qm0pW2ajBnbrWlVy1TVPSJyD65yyzIRKfgnTMcCT/hhu4DjgF0iUgCWAs/VONbNwM0Aa9eunZMbvUYxMtRFYekwe4Z7GNrfQ7z4wNQ7zQH/fup1nPL1K9CiMrZMEQ2Iu4TCQSAU9h8TsvA3ERrA2JICqKLFgP3HdXPvZxtbaOSxD6xn9XUDBCVBY4i7fChl2XnfFIi71T1JGRMoKFIOiLtjCNStE9IKmoovoBJnwg78Z3+tSuzEYxLm+di5s6v+2dK051OmibhfRP4CCEVkDfAB4N/m8gsaZZu6+47TuZGhhtFBtJl9UtUviMhPgNfg/ud5s6o+OI/f3xD7tESWm3UyjDxtZpsmYzZ2q55qmUf6p074Znp/DDwIfB+nIgHOAL7h39/pP+O3/7OqtrUBklAplQrsf3oRS5ce5O5Xz1+u1wNvvoLB956PHFZibJEiZddfDoHyAij3BkTdAeWegKFjeyj3hA0Xdglxl1I4KL4XHYTDQtyjRD2u4EncGxMviIm7XBGUaFFEMObEH+BEnA/HDHyQigvXzIqmaKh+/+yBjERCUOqopzPjmWUolYi8VkQe9pXXNtTYPmFlNhG52K9/WET+dA7O5v3AC3EFBr6ICzf64GwPOi+2SWl+wQtbbGnkMhPmMsyzwYhIICL3q+pDqvpJVf3EfAg7u3cyjCbQRrZpMmZrt+rx3K0AbvdxnwHwZVX9pog8ANwhIlcBPwNu9eNvBT4vIjtwT51On9YZtSC6u4uxqJsj1zzL07uW0XfzZnaefcG8zuGxv7gEgNVf/hjxfy1wAmeUeRNytXj0vPX0f2IL2qV0PRO6RuYRxD1OhAYjARqq64VXJBN1AiK4iy1QFMm86T5HL8nJQxRC0LJAUSEWH8MprBrYWtFgvWNQZhUK6K/VTwJ/gnsafK+I3KmqD+SGpZXZROR0XGW2PxeRU3DX7AuB5wPfE5ETfXjRTHmDql4KXJqb41uBf5jFMWG+bFNnhXkYxuyYpX2ab1Q19u0IjlfVX8/jVx/y906GMa+0mW2ajNnarSnFnar+AnhZjfWP4ookVK8fAd463Ym0Ml3PhZQXxew/2MOy5+3nP/7sY5zx4zO5/eW3zftcdrztskm3r9q2lUfPmz/B07U7RGIYPTKisC9ASuIEmPe+hcNBWgVTQyUYDVyIZkR605y0Tkhz8nz9FJcAK1Am65UXqvMUFhUVmffznTdm95Tp5cAOf40iInfgKrHlxd1puEIn4CqzfcJXZjsNuENVR4HH/I3Gy4F/n8V8Lma8kKu1blrMl21q00Rsw2gcbfYUHCe0fikiPwaGkpWq+qZGfaHdOxlGE2g/2zQZM7Zb08q5O1T51WXn8cIN2ziwtJuHz7wcgL2l1qzmMd9Cp2sv3L/lPFZvHCAsCeVeJShDXBBQpbwkcoJPcI3Pg1whFS/yXIXNXK6dL5qCum0SuUXJRCMRIEpchP6/3cLgezsr/06mNlBHiMh9uc83+1wMyFVd8+wCXlG1f0VlNhFJKrMdA/ywat9jmAEi8jrg9cAxInJDbtMSoFx7rxajDcM5DKPR1GGfWo2PNnsChmE0nja0TZMxY7tl4q5OSkuV4sIx+j+3kcF3buDrf/i3zZ5SS3D/lvNYc80AOy5Zx4lXDbhwSoGgDKIB5WKUCrVoYUww7EMzk951SXimaJoBKomnLiDtb5eUuFXJqmeikorDjhJ4ivN+Ts4zqrp2gm31VF2baExdFdvq5AngPuBNwE9y6/cDbdOB3jx3hpGjPvvUUqjqD5o9B8MwGkyDbZOIvAS4EVgEDALvUNV9InI4LgLq94DPqur7Jtj/rbiIqZOBl6vqfbXGJczGbpm4q5PSEWWWLxpm7xPLOekfr2Rkdw87z7qw2dNqCYr7c7l0ZYh7XYHBNJRSXRXNcMjFW0a9vn/dqCvEkog2vNjTgMyLFyiSVND0YZva5StsBjGEggbKznPmNwey4czu6VNSdS0hX5Gtekx1ZbZ69q0LVf0P4D9E5Is4W3O8qj48k2M1Fcu5M4xK2uzpuIjsZ/ys9+IePq1PQtgNw2hzGmubbgHOV9UfiMiZwAXAZbgm45cBL/LLRNwP/E/gpnq+bDZ2y8RdnXQ9XWDP0HJkxQjBzxZD/xgnf+0KDj69EGDeC6y0Eg9c67x30QnDBIO9xN2KBoqoEA4FTsApxAWcWCtL2tYA8a85Dx2SVc1EnSAUzYlFfDsE9a9CU4rcNJJZeovuBdaIyErgv3CJ+X9RNSapzPbv5CqzicidwBdFZABXUGUN8ONZzQZeC2wBuoCVIvJS4MpG5rvMGWqeO8Oopg2viQHcQ6ov4v4XOR14HvAwcBvwqqbNzDCMOaPBtukk4F/8+7uA7wCXqeoQ8K8isnqynZNql1U9fydjxnZrylYIhmP7JetcWf9nuxl+fkS4P+Tg0wvZefYFqahY+cVreOGGxjY3b0X6btzM9kvWwX/1EhdxXrmyuBBMX0o+7lEINK2Q6ape5g4SJGKO7DmFyvj46aprQmL/HaHS/6ktDTrDJjCLVgi+f9L7cIbnQVyVtl+KyJUikgiqW4HDfcGUdcAGv+8vgS/jiq98Gzh3lpUywYUhvBzY47/j50D/LI85f0z1W9hiSzsvM2Guj9d4XquqN6nqflXd5/OTX6+qXwIOa/bkDMOYIxprm+7HpZmAK3503CRj54IZ2y3z3E2DcEkJfaoHDZRwWCgcN0zfjZvZec4F9P/dtfRs7+WXG8/jBZdv46Er2yalaNbsPOcC+j+5BSm48MlgJHAtEbr81SQQHHTPEbSgzltXUCfywF94mWoTBY19u4R0pR8Tq2upUMT3aRLoUigLg3/dGTl3orMv56uq3wK+VbXu8tz7CSuzqerVwNWzmkAlZVXdO42nVS2D785hGIZnLuxTE4hF5G24vBjI+sxBK0tSwzDqpg7bNFkhOncMke/hvGPVXAqcCdwgIpfjop9Ks5zyVMzYbpm4mwbBzl7KS2KKewJKyyN4tpfCshJ9N26meFhE1KOcfOk2xl58kL7briPcXUSOGeaR0y+tOM7L/ulSdj+7mMF3just3bYMnns+/Z90njONhXiB942HSnAgTPPlgpJ4YQeuQqZ/H2vaDF29Z0+qnra4RuduTIoPy0x66HUKbRj2NBn3i8hfAKGIrAE+APxbk+dUNx32WxjGrGnDa+IdwPXA3+L+V/kh8Je+uXjN4geGYbQfU9imyQrRAaCqfzzFV5wKICInAm+Y1uSmz4ztlom7abDjonWs3jRAeaESjAQEZUF29xJ2KWN0010Whl8wgg4VCQ6EPLp+Hf/tW5em3r2En71hLp0i0HfLpnkv7rJq21YKBwK0oIwtjhl83/muaMqw897FUeCLpgRot8KYOE9bwXne4oK6huWFJAfP9a5DSasdqY/cTHLvtOAEnoaKiiKBpO52DTrs4Wtnnc77cU+9RnGx498BPtbUGdVLa4eaGUZzaLNrwhce+LMJNv/rfM7FMIwG0kDbJCJHqepTIhIAH8ZVzmwYs7FblnM3TXZcuI7HPrDeFwhRykti4h4n9uIuRQIl2B+iRaXvps3s+eXhEMApX7+Cvts3svL6rfTdtHnccU+6Ynq5en03b049ZZSFvhuzY678wjX03Zx97rtlEwD9f7N1Bmdcm0fPW8/o88eIut2VdMKWAbSglBfFaFGJe2MIIO6NCA4GTnwJxEV1+wQKoXOha0HTwihJMRVCd1wtes9eQdFiTOy/T1TS6pptGCI0OeqePk22tBmn+KUA9OAapd/b1BlNg6l+C1tsaedl2kxhn1oRETlRRO4Wkfv9598RkQ83e16GYcwhjbdNbxeRXwEP4QqdfCbZICKDuAIo7xKRXSJyil9/i4is9e//h4jsAl4J/JOIfGeyL5uN3TLP3QxRAe1WCvsD4qKiR48yNhrCUJGug0JpQUy4PyTuUsIlJQ4+uQgKMY99cD39n9jCC7/xEX55Wtaf8OEr6s/RW7VtK9ItSCSsvH4rEgRIDKsGthItiUAL0B3Rd9t1IHD0sXsAkOWjc/pvsPM9F/KiOy/nwONLiBbEpP3sFkVpDp2MBk7ouZ7mBL6heVzEed8SL13oc+nEvXc7ayrgEtJtgCDuGL6iZkfRWefzBeB8XDJyi97+TYK1QjCMStrPPn0aV7b8JgBV/YVv0XJVU2dlGMbc0kDbpKrX48Ika23rn2D9Wbn3XwO+No2vnLHdMnE3Qx770HpOunIbxKALQXb1kLRiiwvuTXFICEYFnlpA1AOjz4/pu2kzQXlmDtNXfGcDTz14pPNsxSBj3s2V9IYrKDISuBBGL4KkEPPbJ5bR9+lNSHdA/99dy+BfXjxn/w73v+lKJyoPK4NA+FyBqCA+HDMmHHLiVyIh7o6Jxcda5vLkUi9cQOVr0iJBSfdJ3yf/BoE/RrmzbsA7rIjH06r6f5o9iRmhreuNMIxm0Sj7JCLLgS/hqukOAm9T1d01xkXAf/qPv66jrcoCVf1xVVGn8qwnbBhGS9Fh904ztlsm7mbBw5c7b9uJVw2kTbc1gLAkyHMFJIKRI2Pi3tjlo+0psHDNHoZ2LGX40SXT/r4f/elGVj64lbhHkZKkf8VBSZz3EHHhjBEQiRN5eSE5VEALc3unumpgK9HhY+w8Y4MLDU28aD6XLk7+wgRkLHBzrspjUtGs70favNyN0TDz7qXhl+r2SUUeMPi+zqiUmdJZBuojInILcDcu7w4AVf1q86ZUPybuDKOKxtmnDcDdqrpRRDb4zxfVGDesqi+dxnGfEZET8DMXkbcAT856toZhtBadde80Y7tl4m4O+NWH1/Hi9duIi1BakgmuqEspHBSiKEQUunYL+55YDEsjCs8W6P/ElrpFyaqtA64C5eFjyO6u1FsnkbioMQXx0ZCiQDlIwx614MMbI2COy9E/um494PP5un2OXDF2otIXOZHRAC36IihlJzoJnDcvFW9JVwTf/DzNp9PcmCREMzkF9d7KWFh5w1Ye+8D6OT23ptF53qJ3Ay8ga2AB7hdvC3HXYf9ZGMbsaKx9Oo2sMe/twD3UFnfT5VzgZuAFIvJfwGO4SnSGYXQKnXfvNGO7ZeJujvjPrc6Ld/Kl29AAyouVYFQoL3Reu2hBzPDzhMK+kKhXkf4hijsWsvKGrUgExX0BD39k4ry7R9evo+/Gzcj+otNGSVtpyVoBRIsjZMy3F4hw3jBx4sp5zASJof+TWxg8d3qerr6bN6fN2mviq1/GvZEvhhK719HQeeEiEA0Ixly+nWt7oNmNcyxoMXZiTlyxlQSJJauOmVsHWTimdFCAjdBxBuolqvriZk9iJkjn/WdhGLOiwfbpaFV9EkBVnxSRoyYY1+P7VZWBjar69VqDRGRd7uO3gO/j/kcaAv4XrgCCYRgdQKfcO82F3TJxN8c8eLUTaC+4fBsohAXnWevaE1BaHhN34So/Di6kMCL0PCupMEyYSEj1/KbAyIpy6q2TpPm3bylAoCBS0fPN9ZHzIYw4MTVdYQdMKOxO+scriR9eRGH1MNETC4gXe/dhqFAKoDtCR0JXGbMMsS+Hmfaq85447Yorwi+T1gZpsRTxoZj5dUluHk4wdhSd5S36oYicoqoPNHsiM6KzfgvDmD2TXxOTNgqeoklwvRyvqk+IyCrgn0XkP1X1kRrjFvvXk4DfA76B+1/nfwP/Mo3vMwyjHeiM/69nbbdM3DWIh650gu2kK7elBUXCA65dQvfTISPHleja18X+/phTLt7GA9eex4vXbyPqhp3XXMBL3r+N//ibTPT13XYdwbJCJtx8dCOiLhQSkOEQiagsTJLGbLoXFk/PxdX/uY3owbCiT1/Cqq0DRMt6KKwapryvm57jhhj9zQJkWQnd14V2xTAcuvy7kcDNuRhD4L2Jvv2B+jy8ijw6fw5JNUyJhXwhlTQ3zwvCwfd3SEgmdGJowR8CZ4jIY7icOwFUVX+nudOqjw5L0DaM2TG1fZq0UfBkTYJF5LcissJ77VYAT01wjCf866Micg/wMmCcuFPVj/rjfhf4XVXd7z9fAfzDpGdhGEZ70SH3TnNht+ou2ygioYj8TES+6T+vFJEfich2EfmSiHT59d3+8w6/vX9aZ9VhBCUIh50AKRwUuvYFEEPvYBfhCHTvDtAAXnT+NoJRKC1Tl7/mbyhXfvEa+m7fCJEQlMQ1/E4KkqgPG0tCHJPKmYk3LCmm6fPVtCtm8J0bpjV/jQRZVK7om5fw6Pp1MCZEz3QjXTEju3sIlpeIh4puPpFk+XShn2/kPYsFtz1tZRB4jxyZZ07UjXVNzJ0X0nkiyfLuCtoRF/M4dIqlvXgtsAY4FdeQ841M3Jhz2jTcNk31W9hiSzsvM2Guj5dxJ3CGf38G7ol1BSJymIh0+/dHAH8ATBUVcDxQyn0u4SpyNhy7dzKMeaRxtqkZzNhuTcdz90HgQSAp83gdsE1V7xCRG4H3AJ/yr7tVdbWInO7H/fk0vqejePCqzPv2go9so9ytaWXNqBsnUmIhGIOoB6JuJ1b2nhyx8gvXEDzRQ6gwdliZ8tElZH+BYNQVK9EApCRES12em5QCJ/4ALSrhqBB3afpHHe6b+Ode+YVreOwdl4zfMBYQ7CkQe+9g36c3QSTsPOcC+m7aTLishEZCPFxAuiMKj/QQrSg7VVl27RAY80VRkLSZeSrwBBR14jPUNMyyov+d+txBLxSTkE3E59y150U7KZ0kWFV1Z4O/onG2STvrtzCMuaCB18RG4Msi8h7g18BbAXwT4HN8z6iTgZtEJInh2FhHyPfngR+LyNdw/2P8D1zBlvnA7p0MY57osP+vZ2y36vLcicixwBuAW/xnAf4I+IofcjvwZv/+tNyXfwV4jcgcl2hsUx766HnsuGid86b5gigaQNyllJYpI0cocXdMMOoqXMZDRcpLIsoLY5Y8UKTnsW60N3aVJwP3RywAxZjC7oIXTl7pKGmvOzdWnKetBv1/d21NYXfix7ZR2BtS3BegiyJWXzfAzr+6kHAoZNW2rdATER0oEg8VXR7hSMjokZHz2BVjKMauJQNAqK4Ng2RzS6tg5v86fBNzDTTNsUt63SWeu9TDNyYMnns+j32og0IyoXFP2zuQebFNzfas2GJLI5fpMtfHyx9a9VlVfY2qrvGvz/n19yXNgFX131T1xar6Ev96ax3HvRpXtXc3sAd4t6peO7vZTo3dOxnGPNJA29QMZmO36vXcfRy4kCzJ73Bgj6omCVy7gGP8+2OAx/3EyiKy149/ps7v6nh2XLiOEzYPEJYExmB0eUxxnxAXITyYNMwj9XD1Pv8A+1kEAQRDoROAIy5/TyJBhgpEPUo44vZVceIxWhATHgyIl5fhQDjhfPRggb6bN9P1TEjpyDLFJSV2/PmH+dVl59H/N1uJFsVQFsrHjLLq768hOjx0YaJ7inDkKCjEQ0XXVF1BuxQd8+fh8+pSkVeISXvh+VYIxL7lgfrWDWHuKgxwotULvDTnMHTbVl6/lcc+2GHijo57+tRIGm6b7LcwjEra8ZpQ1Z8CP53nr7V7J8OYR9rRNk3GTO3WlJ47EXkj8JSq/iS/utYc6tiWP+7ZInKfiNz39NNP1zXZTuKRC9YRF5SoV+l5KqB0RIREEPe4deHBgGDhGDImxL9YmhYeibvjVDBJJJQXRQRl3DqvmwhcNcpwKCTuVoI9TsP3f/7aity5l/yfy7IJiRKvGkZKATv+/MOs2jpA/+evRbtiL64EHSoQ7+ki2FuAQkxw9AjxWEB8oOjCQSNBe9z8ZNiLSYW0eqYojAXZUxR/DmmrAx92KWOSrpeSX3INzNP8vA5+pik6+WLMj22KDg5N+VvYYks7LzO69ub4eJ3IfNinMUbnYKaG0TmYbXLUE5b5B8CbRGQQuAMXUvBxYJmIJJ6/Y4En/PtdwHEAfvtS4Lnqg6rqzaq6VlXXHnnkkbM6iXZl+yUuTDJa4PLltKgEwwGiTuTpc91oqIweGaG9kcu1AyhnRUYKB0KIhcKBgKjXtUSIC7jX3pigJGi3UjgYwJ4uCJS+WzcBsPfxpay5doBlK/YhZe/1WxjRd9t16PNG0KEC4VAA3e5RiCwsQ4TzrA2HlPd2OY9c4L93YeT+e/INymU4RMYC974UkDQvd1+UePni1HOXeufKvnDMmHtNKmcGpSRZ0Y0jdufZcfhzm3QxYB5sU9i7cOrfwhZb2nmZLlPZJyOh4fapSHdjz8Aw2gmzTSlTijtVvVhVj1XVfuB04J9V9R24pnpv8cPOIKtqla929RY/vgPvwOeGHReuIw6heEAIxlwBFClT8YeogSIjIfGyMZK+cG69GxeOkB5DItDuGBkTglEhLoCMuvfaHRMcCCk+U6D/E1vQ3ojS4RF7nlzivIK/7SHYW0BGQnRvF/RGxL2KDDkvnJYCF3JZUGRBlDUQLwVIyYuzssBYQLy4jPbEaCFOc+goxM6DFyV96tSN903O07BLpcK7R+T+bSQiC8v04/I9/ToF3+HCnkBNwXzYpqQpqi22dOoyXaayT4bD7p0MY34x25Qxmz53FwF3iMhVwM+AJKn5VuDzIrID99Tp9NlNsfN55MJ1rLl6AHyunMTi/uONAMS1G1gQQykgOBiQCLxwRCgvjtERJ95KK8Yo/raIr9VC3KUEIwFBGcoLYwp7Q6KjS4yNhNDljkfge9AtKaOxE4AyGqA9EYyGaG8EIz6UMlQoAz0xDBVc4/Hh0FXEBFdAZTSsEGWAC9fsjlxIZkGdJ7AskLRJSFr3Re68k3UaeAErSlxwLSA0GZu8n8HNSTtwqBmiOWZObZP9FoZRiV0Ts8LunQyjQZhtckxL3KnqPcA9/v2jwMtrjBnBly826mf7petYc80AEuPaDigkDbqjpWWCA6HLNYtBC07YRd1KOBQQ9cZoCDIUEvUohX1heoxocUSkXjCOQvBMkWhpRNhTRvf0EC+IXZ6cQjgcEPVGTtiB86z5UEiJBC35xuMjzlOnRWBxGRkqoIUYCTTtXcdokIo2Db0S647T4wUHQuJlZWR/wQta/w/hhWtcAEIlGBG04NsiJE46H8IpKgTlzvPcAR0rWhtFw2xTEuZhGEaGXRPTwu6dDGOeMNsEzM5zZ8wx2y9ZxwlbBkAUKYsTaALhvhCJhbjgWyCUnQAMSuKqZA65ypl4T1d5UZQ2C9fYFWeJFsQIQnT4GHKwQDzcnTYWpywUdndRPnyMwjNFykeMudDMpCedis+Nc8VNVBTtciGVGgp0R0jBCbi00qXiPXmZyCMRh2NC3OO3BepEmw+z1KK6XLsYfzznwVSBtCh04mJPvHydxiEYQtDKiEVGGUaG2SfDMFoRs00pJu5ajEfOX8fqTQPEBVeUBNE0TDMoi6tDkoi8UAnGvCCKST19lF0lyXhhhBwMiXs1LUYSPlsk7lUIQUMXgkmgRCtGYTSkvDRChkJ0YeSOFye5fcl7d+Uo6kIsjJFFGgAAIABJREFUSy5kU0u43Dnx65N5dMVO0CXFVMZ8gRbB7es9dYoTkgRZkZQkn04lJ+Z80/c0LrNDHXczyYUxGoP9FoZRiV0ThmG0ImabHCbuWpAdF7oqmidsGkALpKGawag471XkvXpe5EhZnADDjdHQhW7qmKShjeGwK64SLYh95SDncVPx3rm9RcKSEC0ro6EgBwqpSNPAedMIyDxlEWnOnC6MnFAr+lw63BwpqGtmvr/oRGFX7MZ4oeiyXH0YZ0Aq9AQvJLsUIgjGkkQ7/92+NYPEdKy460iPZDui9p+FYYzD7JNhGK2I2SagvlYIRpN45MJ1TtgorqVBgPPUQXrTKV7oifrwRd8DLyg5r5aGTvyVFyjlI8bQogvTTHvllQNkJEB7I1ep82AIY97LFuFy27x3DkibjqdiLMDtk1TBjN0iXU7o6WiILii77fh9kovPz48gdzxIBSdJaGYi4OKckkt87514ISe/7RxWuDNmgdpiSwcv00XNNhmG0YKYbUoxz12Ls+Oiday5diANTdTQ/5FGQmHYiZ2xAOdd89UmEzEUjASuWEmy33BIeDBwrQmKMcG+gguBDN3xJBLU5/kBmRBLQiAjdf3o/CMBLcRuP5wYlFDR2FX31FiyurTlIAvVVLL3AKruBJKv9CGlyY1Hdfy05p144sJYO42k/L7RfOy3MIxK7JowDKMVMduUYZ67NmD7xevYsWGda40QuzDFcMQVWIlDJUh6HyiufUBBiX0/urjoeuCJAuLCMqOlETIWOGHne+UlXjENsobqiWeuAp93p11eRSbet7KgidevqFnFokS4jWVNyp1YxHn6fNNzAh9emvTES3Lxgpy6S96L++q0eEsnMtdP23OIyHIRuUtEtvvXwyYYd4Yfs11EzvDrFojIP4nIQyLySxHZOLvZtDjN9qrYYkujl5kw18czDMOYC8w2Aea5ayuSdglxUdEQgjJOzEHaEkAVgsjl3ZW7Y1dls8cpLQ195cyCuEbko4Hfgcxrpz4V7qBrWI6SeuAkSkIk1Ym1Ll8YZSwA9aJO8c3MvahLvHBJ2GZREQEdy+Xo+eOnhVR8KGrqrvRNzjURmv6Yj563vvH/6M1AQeKGWqINwN2qulFENvjPF+UHiMhy4CPAWjcjfiIidwKjwBZV/b6IdAF3i8jrVPX/NnLCzcSeBBpGjsbbJ8MwjOljtinFxF2bsf2SdZx05bbUoxV3QRCBhrj8OHDFVMSFZWpSdTPJz+tRKMYuTy5QJ+hiEO+5kzHvtYOsOqV/n4ZERk4cuvYFkj0R8dU9K6tZJsfy+Xj40E3vrcsqfAaZkExEnLq5OW8eoK6R+WMf6FBRl6PB5XxPA17l39+O6790UdWYPwXuUtXnAETkLuC1qvr3wPcBVLUkIj8Fjm3obJuMiTvDqMTKjRuG0YqYbXKYuGtDHr78PABWbxxwws73vkvz8cq4X9bn0RG4MM0k7FGGQ98KwfeUK4tvr+C8aqn3zYdG5ptCioLGyXapHJPrQ1eZHEel904U6YrRsmQCL4kbTdznyXErwkIV7djymJXUISiOEJH7cp9vVtWb6zz80ar6JICqPikiR9UYcwzweO7zLr8um6PIMuDPgOvr/N72Q7E+d4ZRhT3wMAyjFTHb5DBx18bs2OCKiay5eqCyTUFA2iYg3zNOA1wophd1GvpG6AtipBSk3juSIiyBZsf04ZgakhNe3jtXiDPPXNKHLiERauD73bn5aeQLs/jm5ck+bp5+zsnxFAbPPR+A/k9umcN/wRZmaj3xjKqunWijiHwPeF6NTZfWOYNaKjqdlYgUgL8HblDVR+s8Znti2s4wKrFrwjCMVqSBtklEXgLcCCwCBoF3qOo+ETkc+Arwe8BnVfV9E+y/GfdAvAQ8ArxbVfc0Yq4m7jqA7Zc6kbd640DmxStk4YySeMgSQRZ5YTccEHcpwWiQFS5Rcbl8oRNXccHl7SWhm4q6XDzwglJdSGVSJEXJhWaShW2Kuvy6ABijUjrkRByCb96ebMuEHVS+71h09k+fVPWPJ9omIr8VkRXea7cCeKrGsF1koZvgQi/vyX2+Gdiuqh+f3UxbH3sSaBg55sA+GYZhzDmNt023AOer6g9E5EzgAuAyYMS/vsgvE3EXcLGqlkXkOuBixqfEzAkm7jqIxJO3+roB9/BCcK0QksKWSX2UyIm/REipQDgSuMbo6hqmS9kVWMkLLc15AoG0Dx3FuKKdQdaYjiz3zs8n65GnmZcv6Z2XD8s8hJ8Mu3K+Df0HuBM4A9joX79RY8x3gGtylTRPxRkiROQqYClwViMn2Qok3TwMw3DMg30yDMOYNvNgm04C/sW/vwt3n3SZqg4B/yoiqyfbWVW/m/v4Q+AtDZklJu46kh0XZb3fVt6w1XvS/BONWLL2B2kBE+/hi/CVWPzOSeETsnESC1qM3fHGXBEUjYOKJuRSdqGVEovz8kVkqECsLlEw9Ra6fWUsU5F5TXgo0mBBsRH4soi8B/g18FYAEVkLnKOqZ6nqcyLyMeBev8+Vft2xuNDOh4CfigjAJ1T1lobOuFmYl8IwxmEPPAzDaEUabJvuB96EeyD+VuC4WRzrTOBLczGpWpi463CkjAu1THLwfN5bRcPzgm+pEJLeyQaRoImrz6MFl6uHF2EaaubaUN8awRdnkfz+oq5x+VjAOE9e5N/HICUnLNWHjg5+qPOrYtZEvdBu1OFVnwVeU2P9feS8cap6G3Bb1Zhd1M7H61jMS2EYORpsnwzDMGbE1LZpykJ0U9QrOBO4QUQux0VAlWYyTRG5FFf68Asz2b8eTNx1OI+uW88JmwYqqk5KBHF3FnIpJaksvFJUot44y3/zAkzKbv9wOCDu1oo8PYnEefAgK+4Sknn/yvnYTj+RXI87VIgXxATDrmdemtfn6f/kllyLBBh8b4fn3pmeaB3stzCMSuyaMAyjFZncNk1aiA4mr1fgORVARE4E3jCtubn9zgDeCLxGtXGluE3cHQI8cuE6Ttg0QNyd5NRBUHIiS0OtbDfgPXyiAmUv7jyKy4+LeuNMwImvwClZGGbScoGk2mZSPTNUpBQ4j152UBJHkAv59L3tAqX/U1sq8+98vl6ybfCvO1TgqXmLWgYLyzSMSsw+GYbRijTYNonIUar6lIgEwIdxlTOns/9rcQVU/ruqHmzEHBNM3B0iSAzBqLjqmGUfhtnlPWqqqcdOQ5cjl/S+y/ewE3XFMQlcg/R4YYSUnLDTrqxZOjHQpZmoS0Ixy4IW4iwUEy8iIRVwGmiav5eEdSbVOpOKmnIIJONZTkuroNbnzjCqMPtkGEYr0mDb9HYROde//yrwmfR7RQaBJUCXiLwZOFVVHxCRW4AbfdrLJ4Bu4C5fr+CHqnpOIyZal7jzk96PK41RVtW1IrIclwzYj+v38DZV3S1uxtcDrwcOAu9S1Z/O/dSN6ZBU0kw4YfMAsfqilVFA7PPpJMrn37k8u2BUnAMtVzQl7o6dqBNFAie+NKmACVnbg8Qzl4i0XB+7RMBJMjaPjhd++fed3BLBVXxq9izag0bbJrH8IsOowOxT/di9k2HMH422Tap6Pe4arbWtf4L1+ToGk1bTnEuqb6kn49Wq+tJcvOoG4G5VXQPc7T8DvA5Y45ezgU/N1WSNueORC9YRjgpSdoJOIiEoubYI6v8qJEpCN8nlxpFV2RwLCEoBuiByjdMjcQIuxIVgJg3Kq8SZxJIuSVPzBFFJRV0avinZvhVhnp2KKhJPvhgVNNQ2TfVb2GJLOy/TZgr7NBtEZLmI3CUi2/3rYROMO15EvisiD4rIAyLSP6svbix272QY80EDbVO7MR1xV81pwO3+/e3Am3PrP6eOHwLLfKNko8WQGBeiKZl3IigL4agTesRZgRU3vqoPnTgPnhwMK/LhNFSk7BqjjwuhTJxxgff66XgBmBd1FQLPL4dESJBOsRiTMbe2aarfwhZb2nmZCXN9vIyJhE81nwM2q+rJwMuBp2b9zfOH3TsZRqNonG1qK+rNuVPguyKiwE2+dOjRqvokgKo+KSJH+bHHAI/n9t3l1z2ZP6CInI17OsXxxx8/8zMwZsz2iytDNVcNbAWVtMiKRKCx6z+nRc2qZvoqm0lOnnovnQaa9qvT0F9J/qKSfJVMqvLmfC6diuZCOck8dn578vrYBzq8RYKCRIeYJZo5DbVN3T3LsBA0w8jRWPt0GvAq//524B5cAYIUETkFKKjqXQCqeqBRk5kDGmqfeljQ4OkbRhth904p9Yq7P1DVJ7wRuktEHppkbK2YuXH/2t7I3Qywdu1a+zWazJprBwhCHw7pwx7joqaRmGkTdLywk1xsc1o3RVxmQRLWGUvqvVNRxOfeVYRZ5hAVKirDJr35RCGEnedcMMdn3cLYFVEvDbVNi5ceq1hBFcOopHGXxETCJ8+JwB4R+SqwEvgesEFVWzE7tqH2aYksN+NkGHnsigDqFHeq+oR/fUpEvoYLg/itiKzwBngFWVjELiq7th8LPDGHczYagCbes7RSpQvRjH1rgoriJlJjx8S75reL5gqnkAvphFSwCblj5o+bq4qpuSu178bNiErntkDIcajFh8+UhtsmexJoGOOYwj5N2ih4iibB9VAA/j/gZcCvccVJ3gXcWuf+84bdOxnG/GL3To4pxZ2ILAQCVd3v358KXInrzn4GsNG/fsPvcifwPhG5A3gFsDd5Eme0Lkm3gswL51+TnnWaGxODBPgPoKmyo0LASb7wSU685fvcpXl3iYcuv4vvnZfuF2qlZ6+DOSTyCmfJfNkmC8s0jEqmsE+TNgqerEmwiEwkfPLsAn6mqo/6fb4O/D4tJu7s3skw5h+7d3LU47k7Gvia78lQAL6oqt8WkXuBL4vIe3BPz97qx38LV8p3B66c77vnfNbG3JN43GKvz9ICJpVXSl7gEaT6LiMp0ZP35OW/o6qEz7iCK7Xy8vx3Ja99N29m59mdG6Ipak+f6qThtsldAvZbGEZCg+3TRMInz73AYSJypKo+DfwRcF+Ncc3G7p0MYx6xe6eMKcWdfzr2khrrnwVeU2O9AudWrzdam+2XrGPNtQNoAOL6mjs088ppvlJlvh6K1hB5eXLFUJzXLxucD89McvPSxunJdwU5L56SFWvpYCwUcGrmzTbZT2EYFTTQPm2khvARkbXAOap6lqpGInI+cLfvDfcT4NONmtBMsXsnw5h/7N7JUW9BFeMQIGmHoEGuQGVS1CTJxZto53zuHdQemAv7TMMxwRVdiaUyLDNXNbOiHYLvy9f/N1sZfH+HVs08BMv2tiyq9p+FYeRpoH2aRPjcB+SbAd8F/E5jZmEYRlti904pJu6MFJcjpxWFVQBQJ6hSj9lEXjrNufzyY/KFUnK5dZoTg/niK+nhAtdiQXDCT0WRAN93Tzlh0wCPXLhu3H7tz6HXcLOlsbBMw8hh9skwjFbEbFOCibtDnJOu2Ab4cMvQr1TvsQsqvXdJa4JU/OUcbVkunkANoVYt7NJQzWqhmBwoUIiFwfeeT9+NmzNPXjIPIFocs/KGrYTDQtStPPah6XvyTr50G6Wl2noi0QRFa2Ax/IYxHrNPhmG0ImabABN3hzQnXjWAiFSItQqSQiZpcZR80p3H7zuuQlG1eJtM2OV3DnAXp3ev939iC3TlDhtkx1BRtEuRsRCJhBM2DSDAjimEWv/fbkFKwmMfWs/B/jF2/tWFk46fd6z8fkth4s4wcph9MgyjFTHblGLi7lBGpbI4Sb4qimaaS/NVLvPXTdI2wXvSJBGJUJEzVxGWmWxLdkjW54WeAOUAumK0GxhLuqbnxuUEYrTYJQtKWdBiTN9t17HzzIsA6LtlE4SKjIQM/vX59H32OuiFwfe67RRa1BC06LQOSUzcGUYldkkYhtGKmG0CTNwdkqy5egBwVTETVIBE6FW1L0iqZzrx5kMmq1scgGuloKBxvtwmmShLBFl+fXUIZyIwQ4XIf1dBoVy13ffRy/IAFe3OFGnfjZuRBRE7z9pA3+0b0UDpu2UT0qUMvnNDNoUWrbwpsTVXawnUYvgNoxqzT4ZhtCJmmxwm7g5FEoGViKNcm4J8HzrNC7LEe+fFk/iwSc2Hbdb6jtz7tNUBONE2kbCjxvGSPnf54YnITMSjP56ECoUIROn/u2tdw/UFZTQKxoVj54VeqyBWobG1sBh+w0gx+2QYRititinDxN0hSFCWyuqU+fKxeYdbPswyX0glT05caT7cMi+6vDDU/HrFFV9JGqUn+Xxx5pmrmFPq9dNKEShVY6pDPmNJjy+B8+z1f/5aBv/3xelrS2KCoiUQi+E3jPGYfTIMoxUx2wSYuOtIVm8amLCoyEkf3eaFk4z3fOUrYiYkve28cMrXPlEhFVOp3qoYUPndFZ67XHN010IBH/9ZNeFqb14uLHMc+VYMiVjMz0cUCRWNhf7PX1vjAC2CAiYoWgcL8zCMDLNPhmG0ImabUkzcdQhrrh1AgxqetRxJdUytFkaS01Ah4/LjvBas8PZVCLsg53nLv0KWn5cUXZFcs/J84/O8py5fnTOoUpPVeXuQFWEht626MieKBOryDEXRcf8IrYXY06fWwbSdYVRg9skwjFbEbJPDxF2nUFXF8oTNA2lenMQgEb5JuRNjiVNLvbdO/H6a7yeXONLS0EenlpLWB2m/8lgysRcnYq2qOEvegRZLrqVBlQcvGVvtdasWfrU8hWn1TXdsqRZxiWMv/d5WRc1b1CqoWoK2YVRg9skwjFbEbFOCibsO4MSrBtLWA5JrNF6hmXzOW7VmAio8Yangk9r6xwlErfQQxiABaJRPgKs6fq1qmbW8cOm2xAPoxyTCMW2JoOPHqbh/AL9O82LRi0gRiMtV3sVWQ7G48VbCfgvDyDD7ZBhGK2K2KcXEXQegYf6DW8Zpn7y3Ldmu4rf5wUFufC6qMfHoTaiHctGVbgpVanGyaprJd8c1FGZeK9b01CUb8+cvKJp6KSs2JcVVxAnRnWe0XqXMhEYW8RCR5cCXgH5gEHibqu6uMe4M4MP+41WqenvV9juBVar6ooZNttkoENmTQMPIY0WGDMNoRcw2OUzctREnXjWARMLDHzkvXXfSldvSXzGVM7kQxsSTJ/lCKQFoDEGchWhqlUBLj5OsrziO94zlxybFWIJsbIUQTLx3kBVVgUpRB5kYTCpnJpOrfhqTqNCkPUNabRMIvMDLC8KqfD0Zr/1ah8YLig3A3aq6UUQ2+M8X5Qd4AfgRYK2f0U9E5M5EBIrI/wQONHKSrYGFeRhGBfbAwzCMVsRsU4qJuzZg9XUDTmCFThid+LFtQE7TxLlCKjEVYq3Ca5YTXUkNlGoBphUirMrplu+DlwjBxEMY5vYLspS5tA8eOVFX8Z05d2B+u1ZtS0RgWmAll9OXzxNUP09cDz9Jzyf/j9Lq5BIaG8NpwKv8+9uBe6gSd8CfAnep6nMAInIX8Frg70VkEbAOOBv4ciMn2gpYE3PDyNNw+2QYhjEDzDYlmLhrcVZfN5B5vaJsfVosBR8GqVpZtKSqRkm2IRemGeukoZbjVudEH5B51iTbnqa5pRU083PKFVKpNcF8/lz+RCZrbl4hXqvmo+LckwnJd8dtcP031lt0tKo+CaCqT4rIUTXGHAM8nvu8y68D+BiwFTjYyEm2BBbDbxjjMW+2YRitiNkmoE5xJyLLgFuAF+Fud84EHqZG3o6ICHA98Hrczd+7VPWncz7zQ4QdF7l+das3DlSUZE/EU6KRgkiIQ60slqKCxOq8annPVj4fL87y7jQn0iqEomS5d25F3tuWeec0EW3J/oFS4bXLeQVTAoUoO1ba2JzcsfLNzbVKwJHbJ/HqKRBWeQH9eo3c2FbOt3PexykFxREicl/u882qenPyQUS+Bzyvxn6X1jmLWpJfReSlwGpVPU9E+us8VsNovG2ysEzDqKA++2Rg906GMa802DaJyEuAG4FFuGv3Haq6T0QOB74C/B7wWVV93wT7fwwXORUDT+Gu8ScaMdd6PXfXA99W1beISBewALiE2nk7rwPW+OUVwKf8qzELdmwY35T8xI9tyxVHgSARORWhim5FtXcuTUPzOXepRy8nANOwSrJ1LqfOHSwN66zeH//9eWGXHCfpe5cO1EoBVk8unEzwmsxT3PEl0IpwTvHKt9V73DlBEU016BlVXTvhEVT/eKJtIvJbEVnhvXYrcEamml1koZsAx+LCN18J/DcRGcTZj6NE5B5VfRXNofG2yW5kDSNHXfbJcNi9k2HMGw23TbcA56vqD0TkTOAC4DJgxL++yC8TsVlVLwMQkQ8AlwPnNGKiU4o7EVkC/P/AuwBUtQSURGSivJ3TgM+pqgI/FJFlyY3knM/+EOdXl51Xc/0LPrINorx+kkxU5ValBVLSldk6FSrz+HyD9My7N8ENb86zNi7HLu+9S26YpxJayeZQM4GZX5+fRlVYpmpWOEUjQQpusAQtfrPe+CfjdwJnABv96zdqjPkOcI2IHOY/nwpc7HPwPgXgPXffbJawmxfbpJjnzjDymOeuLuzeyTDmmcbbppOAf/Hv78LdJ12mqkPAv4rI6kmnp7ov93EhDSwCUY/nbhXwNPAZ75L8CfBBJs7bmShXp8JAicjZuIIMHH/88bM5B6OKhz46XvStuWYg+5CIulylSok1C68kV6SlKpSyonVCnsTLlxd2+X3z4ZoVx6zKzav27KXjNfNGprvmiq3khWSUiDzJQkXBt0KggZfTHNJYQbER+LKIvAf4NfBWABFZC5yjqmep6nM+hOBev8+VSXGVFqLhtqknXASReSkMowJ74FEPjbdPLGjoCRhG29FY23Q/8CbcA/G3AsdN9wAicjXwTmAv8Oo5nV2OesRdAfhd4P2q+iMRuR4XRjARNXN1xq1w+UE3A6xdu7Ydbrfbmu2X1A7rzIdS5huTp8KuysGWVtNMtFStgiwxiEjmuQtywivvfav1q6fet3xhFCqFYSL+kmPlvYgVE61y8/mxg395cY0vbiFUGyooVPVZ4DU11t8HnJX7fBtw2yTHGWTyEIRG03DbtLR4lH+QbhgG0HD71EE03D4tkeVmnAwjYWrbNGmtApiyXsGZwA0icjkuAqo0/SnqpcClInIx8D5cy6k5px5xtwvYpao/8p+/gjNQE+Xt7KJSzR4LNCRh0JgdE4V1rt44ULmiOjwTJs6Nq16fL4qS760w1X9JVSIyDS1NRZ5mhVSSMUmFTnL7jGtyPsX3tgomKOphHmyT3cgaxjjMPtWD3TsZxnwzuW2atFaB233iegWeUwFE5ETgDdObXAVfBP6JZok7Vf2NiDwuIiep6sO4J/4P+KVW3s6dwPtE5A5cMvBeixlvL2oVbzlhkxN8tToX5HuIJ8IvDctUcTHQIePdfMn7QF1OXSROsFW3RMh78CqKsTDem1eVW5iOSQRqq9dSAZygsLCnqZg322Q3soaRw+xTPdi9k2HMN421TSJylKo+JSIB8GFc5czp7L9GVbf7j28CHprrOSbUWy3z/cAXfLWnR4F3426Xx+XtAN/ClfLdgSvn++45nbHRFB65cLzgA1j58a3p+4pWCvlcu6B6QBW53L+s/x3jq2iqVDZSz39fru1CRVPz5DpPjhO3gbpTULWbpzpprG1SzHNnGHnMPk0Hu3cyjPmi8bbp7SJyrn//VeAzyQZfQXwJ0CUibwZOVdUHROQW4Eaf9rJRRE7C3ZnupEGVMqFOcaeqPwdquTJr5e0ocG6NsUYH8tiH1o9b1/+pLbUHV/XHy4SbZqIrGZPk6SVCLRIoxFRWd8ntD1lYZkURFkmPr9ri/e3y2JPxumi4bVJF7bcwjErsmqgLu3cyjHmmgbZJVa/HtTepta1/gvX5Ogb/qzEzG0+9njvDqJvBvz5/3Lq+WzdVVrSEykIo1SGb1VU2E+EW5vbJj4N0fwk1i6QTdU3YA239FggJqlaNrpWw38IwMsw+GYbRiphtSjFxZ8wLO99z4bh1fbdscm+SfLmYnKcN39su87yloZl5qqtq4lseKKn3T9X3tmv55uUZaqGArYFVBjSMcTTKPonIcuBLQD8wCLxNVXdXjXk1sC236gXA6ar69YZMyjCMtsHunRwm7oymsfOs8YIPcqIvzoVwQlX+HU685RuX54utjKucSaWnsJVRK1jQKihgrRAMI0dj7dMG4G5V3SgiG/zniyq/Xr8PvBRSMbgD+G6jJmQYRptg904pJu6MlqOW6Ou7dVOuLQLUrJqZ5Nsl5FL0NBZ2vqtN8u0ArGBBi2CeO8MYR+Ps02nAq/z724F7qBJ3VbwF+L+qerBREzIMo42weyfAxJ3RJtQM60zy+KrDLZMCLG2KqlpoQaugFuZhGHnqsE9TNgqehKOT8v++D9xRU4w/HRiYYoxhGIcAdu+UYeLOaFtqCr5Pb6pcUe3haxM0tlDAlsHCMg2jgins06SNgkXke8Dzamy6dDpz8A3AXwx8Zzr7GYbRudi9k0NaIZ9ERPYDDzd7HnPEEcAzzZ7EHNAp5wHNP5c+VT2y3sEi8m3cnCfjGVV97eymZUyFiDwNDNEZ10Kzr4O5pFPOpRXOY67t04xtk4g8DLzKe+1WAPeo6kkTjP0g8EJVPXsm39UJ2L1TS2LnMXe0jG1qN1pF3N032ZO+dqJTzqVTzgM661yM+adT/n465Tygc86lU85jrhCRzcCzuYIqy1W1ZuUtEfkhcLEvsHJI0kl/P51yLnYeRisQNHsChmEYhmEYwEbgT0RkO/An/jMislZEbkkGiUg/cBzwgybM0TAMo6WxnDvDMAzDMJqOqj4LvKbG+vuAs3KfB4Fj5m9mhmEY7UOreO7qraTVDnTKuXTKeUBnnYsx/3TK30+nnAd0zrl0ynkYzaGT/n465VzsPIym0xI5d4ZhGIZhGIZhGMbsaBXPnWEYhmEYhmEYhjELmi7uROS1IvKwiOzw1bFaFhE5TkS+LyIPisipC4N/AAAEA0lEQVQvfSlmRGS5iNwlItv962F+vYjIDf7cfiEiv9vcM6hEREIR+ZmIfNN/XikiP/Ln8SUR6fLru/3nHX57fzPnXY2ILBORr4jIQ/63eWW7/iZG62C2qbl0gn0y22Q0CrNPzaMTbBOYfepkmiruRCQEPgm8DjgFeLuInNLMOU1BGVivqicDvw+c6+e7AbhbVdcAd/vP4M5rjV/OBj41/1OelA8CD+Y+Xwds8+exG3iPX/8eYLeqrga2+XGtxPXAt1X1BcBLcOfUrr+J0QKYbWoJOsE+mW0y5hyzT02nE2wTmH3qXFS1aQvwSuA7uc8X4/rWNHVe05j/N3Dlmh8GVvh1K4CH/fubgLfnxqfjmr0Ax+Iu3D8CvgkIrmFlofq3Ab4DvNK/L/hx0uxz8PNZAjxWPZ92/E1saZ3FbFPT59/29slsky2NWsw+NXXubW+b/HzMPnXw0uywzGOAx3Ofd9Em5Y29e/1lwI+Ao1X1SQD/epQf1srn93HgQiD2nw8H9qhq2X/OzzU9D799rx/fCqwCngY+48MkbhGRhbTnb2K0Dm37d9IBtgk6wz6ZbTIaRdv+rXSAfeoE2wRmnzqaZos7qbGu5ct3isgi4B+BD6nqvsmG1ljX9PMTkTcCT6nqT/KrawzVOrY1mwLwu8CnVPVlwBBZGEEtWvlcjNahLf9O2t02QUfZJ7NNRqNoy7+VdrdPHWSbwOxTR9NscbcLOC73+VjgiSbNpS5EpIgzTl9Q1a/61b8VkRV++wrgKb++Vc/vD4A3icggcAcuvODjwDIRSRrb5+eanoffvhR4bj4nPAm7gF2q+iP/+Ss4g9Vuv4nRWrTd30mH2CboHPtktsloFG33t9Ih9qlTbBOYfepomi3u7gXW+EpDXcDpwJ1NntOEiIgAtwIPqupAbtOdwBn+/Rm4ePJk/Tt9laHfB/Ym7u5moqoXq+qxqtqP+zf/Z1V9B/B94C1+WPV5JOf3Fj++JZ7YqOpvgMdF5CS/6jXAA7TZb2K0HGabmkSn2CezTUYDMfvUBDrFNoHZp46n2Ul/wOuBXwGPAJc2ez5TzPUPcW7oXwA/98vrcTHUdwPb/etyP15wFa0eAf4TWNvsc6hxTq8CvunfrwJ+DOwA/gHo9ut7/OcdfvuqZs+76hxeCtznf5evA4e1829iS2ssZpuav7S7fTLbZEujFrNPTT+ntrZNfn5mnzp0Ef+jGYZhGIZhGIZhGG1Ms8MyDcMwDMMwDMMwjDnAxJ1hGIZhGIZhGEYHYOLOMAzDMAzDMAyjAzBxZxiGYRiGYRiG0QGYuDMMwzAMwzAMw+gATNwZhmEYhmEYhmF0ACbuDMMwDMMwDMMwOgATd4ZhGIZhGIZhGB3A/wP+7y7AkmJpOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 2. read deformation data in meter unit\n",
    "unw_file = date12+'_msk.unw'\n",
    "data, atr = readfile.read(unw_file)\n",
    "data[data==0] = np.nan\n",
    "data *= float(atr['WAVELENGTH']) / (-4. * np.pi) #convert phase in radian to displacement in meter\n",
    "\n",
    "# 3. get lat/lon\n",
    "lat0 = float(atr['Y_FIRST']);  lat_step = float(atr['Y_STEP']);  lat_num = int(atr['LENGTH'])\n",
    "lon0 = float(atr['X_FIRST']);  lon_step = float(atr['X_STEP']);  lon_num = int(atr['WIDTH'])\n",
    "lat1 = lat0 + lat_step * lat_num\n",
    "lon1 = lon0 + lon_step * lon_num\n",
    "lats, lons = np.mgrid[lat0:lat1:lat_num*1j,\n",
    "                      lon0:lon1:lon_num*1j]\n",
    "\n",
    "# 4. dislay\n",
    "display = True\n",
    "fig_file = '{}.png'.format(os.path.splitext(unw_file)[0])\n",
    "if display:\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=[15, 3])\n",
    "    im = ax1.imshow(data);   ax1.set_title('Displacement');  cbar = fig.colorbar(im, ax=ax1);  cbar.set_label('meter')\n",
    "    im = ax2.imshow(lats);   ax2.set_title('Latitude');      cbar = fig.colorbar(im, ax=ax2);  cbar.set_label('degree')\n",
    "    im = ax3.imshow(lons);   ax3.set_title('Longitude');     cbar = fig.colorbar(im, ax=ax3);  cbar.set_label('degree')\n",
    "    ax1.plot(int(atr['REF_X']), int(atr['REF_Y']), 'ks', ms=3) # reference point\n",
    "    #plt.savefig(fig_file, bbox_inches='tight', transparent=True, dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saved to 20160101_20171022_msk.txt\n"
     ]
    }
   ],
   "source": [
    "# 5. save data into text file\n",
    "out_file = os.path.splitext(unw_file)[0]+'.txt'\n",
    "header = 'Displacement data from file {}\\n'.format(unw_file)\n",
    "header += 'InSAR file: {}\\n'.format(unw_file)\n",
    "header += 'Reference point in (Y, X): ({}, {})\\n'.format(atr['REF_Y'], atr['REF_X'])\n",
    "header += 'Reference point in (lat, lon): ({}, {})\\n'.format(atr['REF_LAT'], atr['REF_LON'])\n",
    "header += 'Latitude (deg)\\tLongitude (deg)\\t\\tDisplacement (m)'\n",
    "\n",
    "# skip pixels with Nan value -> irregular grid\n",
    "skip_nan = True\n",
    "if skip_nan:\n",
    "    mask = ~np.isnan(data)\n",
    "else:\n",
    "    mask = np.ones(data.shape, np.bool_)\n",
    "\n",
    "# save to text file\n",
    "np.savetxt(out_file,\n",
    "           np.hstack((lats[mask].reshape(-1, 1),\n",
    "                      lons[mask].reshape(-1, 1),\n",
    "                      data[mask].reshape(-1, 1))),\n",
    "           fmt='%s',\n",
    "           header=header,\n",
    "           delimiter='\\t')\n",
    "print('saved to '+out_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
